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Preface 


The  purpose  of  this  study  was  to  develop  a  variable  wind  fallout 
code  that  could  be  used  to  predict  the  radiation  exposure  and  population 
insult  due  to  a  nuclear  attack  against  the  U.S.  The  operational  codes 
currently  used  are  constant  wind  models,  and  the  code  developed  is  based 
on  Hopkins'  two  step  technique  for  finding  fallout  from  variable  winds. 
The  thesis  was  sponsored  by  the  Air  Force  Operational  Test  Center  at 
Kirtland  AFB,  New  Mexico. 

1  would  like  to  thank  Tom  Hopkins  for  his  assistance  in  giving  me 
his  hotline  locator  code,  the  wind  spectral  coefficients,  and  for 
explaining  how  his  model  works.  Also,  I  want  to  thank  my  thesis  advisor. 
Dr.  Bridgman,  for  his  advice  and  encouragement.  And  finally,  I  thank  my 
wife  Mina,  who  did  the  hard  work  of  running  the  house  and  taking  care  of 
our  family,  while  I  took  eighteen  months  off  for  school. 


John  W  St.  Ledger 
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Abstract 

Hopkins'  variable  wind  fallout  model  is  used  to  predict  the  dose  and 
population  insult  across  the  United  States  from  a  nuclear  attack.  The 
dose  calculation  is  performed  by  two  programs  written  in  Fortran  V  for  a 
CYBER  845  computer.  Hopkins'  hotline  locator  program  was  modified  to 
reduce  its  run  time,  and  it  is  used  to  locate  the  fallout  hotline  as 
trace  particles  are  translated  to  the  ground  in  a  spatially  varying  wind 
field.  The  second  program  analytically  smears  fallout  activity  along  the 
hotline.  To  reduce  run  time  and  to  match  the  population  model,  the  dose 
program  uses  a  computational  grid  of  one  degree  latitude  by  one  degree 
longitude.  A  difference  of  cumulative  normal  functions  gives  the  average 
dose  across  a  grid  cell.  An  analytical  method  was  developed  to  treat 
multiple  bursts  against  an  area  target  as  one  cloud. 

For  the  winds  of  0000  Universal  Time  on  16  January  1982,  a 
hypothetical  attack  against  twenty  five  air  bases  and  six  Minuteman 
missile  fields  results  in  26.9  million  fallout  deaths.  This  calculation 
used  407  seconds  of  computer  time. 


INCORPORATION  OF  HOPKINS'  VARIABLE  WIND  MODEL 


INTO  A  POPULATION-DOSE  FALLOUT  CODE 


I .  Introduction 


Background 

There  are  two  general  kinds  of  fallout  code  currently  in  use: 
discrete  codes  and  smearing  codes.  Discrete  codes,  such  as  the  Defense 
Land  Fallout  Interpretive  Code  (DELFIC),  produce  "fallout  footprints  on 
the  ground  by  numerical  integration,  employing  discrete  cells  in  space, 
time  and  particle  size. "(4:205)  Smearing  codes,  such  as  WSEG-10  or  the 
AFIT  Smear  Model,  produce  fallout  footprints  by  analytic  solutions  to 
equations  which  approximate  the  fallout  cloud  location,  size  and  activity 
arrival  rate  on  the  ground.  While  DELFIC,  as  much  as  possible,  uses 
physics  to  translate  fallout  to  the  ground,  smearing  codes  generally  use 
a  mix  of  physics  and  empirical  curve  fits. 

Hopkins  has  developed  a  two-step  method  to  calculate  fallout 
footprints  which  uses  both  discrete  and  smearing  techniques( 11 ; 12) .  The 
first  step  is  to  locate  the  fallout  hotline.  Discrete  trace  particles 
are  started  at  their  initial  heights  using  an  empirical  curve  fit  to 
DELFIC  stabilized  particle  size  wafer  heights( 11 : 14-15) .  Using 
McDonald-Davies  fall  mechanics(4:212) ,  the  trace  particles  are  translated 
by  variable  winds  as  they  fall  to  the  ground.  Then,  the  curved  hotline 
is  approximated  by  straight  lines  connecting  the  ground  locations  of  the 
trace  particles( 11 : 24) .  The  second  step  is  to  analytically  smear  the 
activity  along  the  hotline  using  a  variable  wind  smearing  equation 


developed  from  the  fundamental  smearing  equation( 11 : 73-82) . 


Hopkins'  two-step  method  has  the  advantages  of  both  the  discrete  and 
the  smearing  techniques.  Discrete  codes,  such  as  DELFIC,  can  use 
variable  wind  fields,  but  DELFIC  would  use  a  prohibitive  amount  of 
computer  time  to  calculate  the  fallout  from  a  nuclear  attack.  Smearing 
codes,  such  as  WSEG-10,  have  the  advantage  of  being  fast  running,  but 
they  use  a  constant  or  near  constant  wind  field,  which  makes  it  difficult 
to  realistically  predict  where  the  fallout  will  land  far  from  a  burst. 

By  using  discrete  techniques  and  actual  variable  winds,  Hopkins'  hotline 
location  method  produces  realistically  curving  hotlines,  and  can 
accurately  predict  fallout  locations  far  from  the  burst.  By  smearing  the 
activity  along  the  hotline,  run  times  are  much  faster  than  DELFIC,  or 
other  disk  tossers,  and  can  approach  the  run  times  for  pure  smearing 
codes . 

Problem  Statement 

To  incorporate  Hopkins'  variable  wind  model  into  a  population-dose 
fallout  code  to  predict  fallout  casualties  due  to  a  nuclear  attack. 

Scope  and  Assumptions 

The  goal  of  the  thesis  was  program  development.  Results  are 
presented  for  a  hypothetical  counterforce  attack  against  the  Minuteman 
missile  fields  and  25  air  bases,  using  the  winds  of  0000  Universal  Time 
on  16  January  1982.  The  programs  were  written  in  Fortran  V  for  a  CYBER 
845  computer. 

Effects  due  to  weather,  such  as  rain,  and  uneven  terrain  were  not 
considered.  Each  nuclear  explosion  was  considered  as  an  independent 
event  having  no  effect  on  the  fallout  distribution  from  any  other 


explosion.  In  the  case  of  multiple  explosions  in  one  area,  such  as  an 
attack  on  a  Minuteman  missile  field,  all  fallout  in  the  area  was  assumed 
to  be  translated  by  the  winds  affecting  the  center  burst. 

Approach  and  Organization 

Hopkins'  hotline  locator  code  was  modified  to  reduce  its  run  time. 
The  modified  code  uses  a  file  of  burst  information  and  a  file  of  wind 
spectral  coefficients  as  input,  and  outputs  hotline  location  to  a  file. 
The  hotline  output  file  can  be  used  to  draw  a  map  showing  hotline  and 
fallout  cloud  locations  or  to  calculate  the  dose  across  the  United 
States.  Using  the  hotline  data  and  Wendel's  US  population  model(21 :5-9) , 
the  dose  program  calculates  the  fallout  deaths  due  to  an  attack  and 
outputs  the  total  number  of  deaths  and  dose  information  for  each  one 
degree  of  latitude  by  one  degree  of  longitude  for  the  United  States.  The 
one  degree  by  one  degree  mesh  was  used  because  it  matches  the  mesh  of 
Wendel's  rural  population  model(21:6),  it  fulfills  Air  Force  Operational 
Test  Center  (AFOTEC)  requirements  for  dose  information,  and  it  uses  less 
computer  time  than  a  smaller  mesh. 

Section  II  discusses  the  hotline  locator  program,  its  modifications, 
and  shows  the  hotlines  resulting  from  a  hypothetical  counterforce  attack. 
Section  III  discusses  the  development  of  the  dose  program  and  the  theory 
used  in  calculating  the  dose.  Section  IV  contains  the  conclusions  and 
recommendations  for  further  study. 

The  appendicies  contain  listings  of  the  computer  codes,  maps  showing 
the  hotlines  and  fallout  cloud  locations,  and  the  derivations  for  the 
algorithms  used  in  the  dose  program. 


Hopkins'  model  for  locating  fallout  hotlines  is  fully  described  in 
reference  12.  His  computer  code  for  calculating  the  hotline  position  was 
used  with  modifications  to  decrease  its  run  time.  The  modified  code  with 
example  input  and  output  files  is  in  Appendix  F. 

Program  Description 

The  hotline  locator  program  works  by  following  trace  particles  as 
they  fall  through  a  1976  US  Standard  Atmosphere (16)  and  are  translated 
horizontally  by  variable  winds.  The  line  segments  on  the  ground 
connecting  the  trace  particle  landing  positions  define  the  fallout 
hotline,  and  the  trace  particle  times  of  arrival  and  wind  shears  are  used 
to  calculate  fallout  cloud  dimensions.  First,  the  trace  particles' 
starting  heights  are  determined  by  an  empirical  curve  fit  to  DELFIC  wafer 
starting  heights( 11: 14-15).  Then  the  atmosphere  from  the  highest 
particle  to  the  ground  is  divided  into  24  (or  twelve)  equal  layers,  and 
each  particle  starting  height  is  adjusted  to  the  nearest  level  in  the 
discretized  atmosphere.  The  north-south/ east-west  wind  components  for 
the  particle's  starting  height  and  location  are  multiplied  by  the  time  of 
fall  through  the  first  atmosphere  layer  to  find  the  particle's  starting 
location  for  the  second  layer.  Then  the  winds  and  time  of  fall  for  each 
lower  layer  are  used  to  translate  the  particle  until  it  lands  on  the 
ground.  After  all  trace  particles  are  on  the  ground,  the  process  is 


repeated  for  the  next  burst. 


Program  Modifications 

Trace  Particle  Selection  .  Any  number  of  trace  particles  can  be 
used  to  define  the  hotline.  The  greater  the  number  of  trace  particles, 
the  more  the  hotline  appears  to  be  a  continuous  curved  line.  But 
doubling  the  number  of  trace  particles  will  approximately  double  the  run 
time. 

The  hotline  program  originally  had  20  trace  particles,  and  this  was 
reduced  to  ten  particles.  The  trace  particle  sizes  currently  used  were 
picked  to  provide  approximately  equal  length  hotline  segments  for  a  one 
megaton  burst.  Twenty  microns  radius  was  picked  as  the  smallest  particle 
for  several  reasons.  First,  when  local  fallout  is  defined  as  the 
radioactivity  that  arrives  within  24  hours(10:388) ,  twenty  microns 
generally  represents  the  smallest  particle  of  interest  for  single  bursts 
in  the  100  kiloton  to  five  megaton  range.  Also,  assuming  a  DELFIC 
default  activity  size  distribution^ : 211) ,  50  to  100  one  megaton  bursts 
must  have  the  same  hotline  to  cause  a  dose  to  infinity  of  100  rads  tissue 
near  the  twenty  micron  particle  location.  And  finally,  if  the  wind  data 
used  is  considered  as  typical,  then  one  megaton  bursts  in  the  eastern 
half  of  the  US  will  deposit  twenty  micron  and  smaller  particles  in  the 
Atlantic  Ocean. 

Atmosphere  Layers  .  The  modified  hotline  program  uses  twelve 
atmosphere  layers,  instead  of  the  original  24  layers.  Hopkins  showed 
that  using  more  than  24  layers  in  the  atmosphere  would  not  significantly 
increase  the  hotline  accuracy,  and  that  going  from  24  to  twelve  layers 
resulted  in  about  a  four  percent  change  in  hotline  location(12:32) . 

Figure  1  shows  two  hotlines  for  a  one  megaton  burst,  one  generated  by 
using  24  atmosphere  layers  and  the  other  from  using  twelve  layers.  As 


Figure 


can  be  readily  seen,  there  is  a  difference  in  the  hotlines.  However, 
this  difference  is  not  judged  to  be  significant.  The  current  hotline 
locator  program  does  not  allow  the  wind  field  to  change  with  time, 
although  a  realistic  wind  field  can  change  significantly  over  24  hours. 
Therefore,  the  error  introduced  by  using  twelve  atmosphere  layers  will 
probably  be  less  than  the  error  caused  by  using  wind  data  from  a  fixed 
time.  If  climatological  winds  are  used,  such  as  the  most  probable  winds 
for  the  month  of  January,  then  using  24  layers  will  produce  more  accurate 
probable  hotlines  at  the  expense  of  doubling  the  run  time. 

Spectral  Coefficient  Handling  .  The  north-south/east-west  wind 
components  are  calculated  from  spectral  wind  coefficients  calculated  by 
the  National  Meteorlogical  Center(12:6).  Using  the  spectral 
coefficients,  it  is  possible  to  calculate  the  wind  at  any  location  on  the 
Earth's  surface  for  any  of  twelve  different  altitudes.  The  original 
hotline  locator  program  read  the  desired  coefficients  from  a  data  file 
each  time  a  wind  component  was  calculated,  and  it  was  discovered  that  the 
majority  of  the  run  time  was  used  accessing  the  spectral  coefficient  data 
file.  The  hotline  program  was  changed  to  read  the  coefficients  into  an 
array  in  working  memory. 

Attack  Scenario 

The  modified  hotline  locator  program  was  used  to  find  the  hotlines 
for  a  hypothetical  counterforce  attack  consisting  of  31  hotlines  for  1125 
one  megaton  bursts.  The  25  air  bases  and  six  Minuteman  fields  attacked 
are  listed  in  Table  I.  It  was  assumed  that  the  wind  translating  the 
center  burst  of  a  missile  field  would  translate  all  of  the  bursts  in  the 


field,  and  that  there  was  no  interaction  between  the  bursts.  All  bursts 


TABLE  I 


Air  Base  and  Missile  Field  Target  Locations 


for  Hypothetical  Counterforce  Attack(9:9;21:38) 


Target 

Number  Bursts 

Location 

1 

Grand  Forks 

1 

47.95N 

97.40W 

2 

Hector 

1 

44 .88N 

93.22W 

3 

Minn-St.  Paul 

1 

46.85N 

92.18W 

4 

Duluth 

1 

56.85N 

92.18W 

5 

KI  Sawyer 

1 

46.35N 

87.40W 

6 

Gen.  Mitchell 

1 

42.95N 

87 . 90W 

7 

Chicago-0' Hare 

1 

41.98N 

87.90W 

8 

Kinchloe 

1 

46.23N 

84.47W 

9 

Selfridge 

1 

42.60N 

82.83W 

10 

Plattsburg 

1 

44 .65N 

73.47W 

11 

Offut 

1 

41.12N 

95.90W 

12 

Lincoln 

1 

40.85N 

96.77W 

13 

Salinas 

1 

33.78N 

97.65W 

14 

Kansas  City 

1 

39.30N 

94.73W 

15 

Forbes 

1 

38.95N 

95.67W 

16 

Whiteman 

1 

38.73N 

93.55W 

17 

Lambert  St.  Louis 

1 

38.75N 

90.37W 

18 

Roswell 

1 

33.30N 

104. 53W 

19 

McConnell 

1 

37.63N 

97.27W 

20 

Natrona 

1 

42.92N 

106. 47W 

21 

Hill 

1 

41.12N 

111.97W 

22 

Salt  Lake  City 

1 

40.78N 

111.97W 

23 

Minot 

1 

48.42N 

101. 35W 

24 

Ellsworth 

l 

4  4.15N 

103. 10W 

25 

Maelstrom 

1 

4  7 . 50N 

111.18W 

26 

Missile  Field 

220 

47.50N 

110.00W 

27 

Missile  Field 

165 

48.50N 

101. OOW 

28 

Missile  Field 

220 

48.50N 

97.00W 

29 

Missile  Field 

165 

44.50N 

107. OOW 

30 

Missile  Field 

165 

41.50N 

105.  OOW 

31 

Missile  Field 

165 

38.50N 

96.50W 

were  assumed  to  occur  simultaneously. 

The  Whiteman  AFB  Minuteman  field  has  150  silos  and  fifteen  control 
centers  in  an  area  of  about  156  kilometers  by  180  kilometers(6) .  This 
field  was  assumed  to  be  representative  of  all  Minuteman  fields.  The  200 
silo  fields  were  assumed  to  have  dimensions  of  180  km  by  208  km,  which 
gives  the  same  silo  density  as  the  Whiteman  field. 


Results 


The  wind  data  used  was  for  0000  Universal  Time  on  16  January  1982. 
The  unmodified  hotline  program  used  1460  seconds  to  find  one  hotline. 
The  modified  version  used  370  seconds  to  find  the  31  hotlines  shown  In 
Figure  2. 

There  are  six  figures  In  Appendix  A  showing  the  fallout  cloud 
locations  in  four  hour  intervals  from  four  to  24  hours  after  the  attack 


COUNTERFORCE  flTTRCK  HOTLINES 


larkers  at  target  locations.  Dotted  lines  show  fallout  hotlines  for  each  target. 


III.  Dose  Program 


The  dose  program  finds  the  population  exposure  and  dose  across  the 
United  States.  The  program  is  based  on  Hopkins'  variable  wind  smearing 
equation(ll: 73-82).  The  first  decision  made  in  program  design  was  to 
pick  the  mesh  size  for  calculating  dose.  It  was  determined  that  a  grid 
size  of  one  degree  latitude  by  one  degree  longitude  would  satisfy  AFOTEC 
requirements  and  simplify  population  insult  calculations  for  the 
population  model  used.  The  relatively  large  mesh  reduced  run  times. 
However,  fallout  cloud  dimensions  at  early  times  are  much  smaller  than 
the  mesh  size,  and  this  causes  the  dose  calculation  to  be  sensitive  to 
minor  changes  in  the  burst  position.  Therefore,  a  method  was  developed 
to  find  the  average  dose  across  a  grid  cell.  An  analytic  method  was  also 
developed  to  calculate  the  dose  from  an  area  target  without  using 
superposition.  WSEG-10  empirical  formulas  were  used  for  cloud  size  and 
growth.  However,  the  WSEG-10  formula  for  cloud  growth  due  to  wind  shear 
had  to  be  modified  to  give  cloud  standard  deviations  on  the  order  of  200 
kilometers  at  late  times,  as  predicted  by  Bauer(18:Sec  6). 


Variable  Wind  Smearing  Equation 


K  * 
if 


£4  - 

A';/- 

with 


burst  point  defining  the  origin.  (Roentgen/hour) 
the  source  normalization  constant 
(Roentgen-kilometers  squared/hour-kiloton) 
the  weapon  fission  yield  (kilotons) 

the  normalized  fraction  of  radioactivity  arriving  on  the 
ground  at  the  time  of  arrival  (per  hour) 
the  time  of  arrival  (hours) 
the  hotline  coordinates 

(jL) 


and 


=  the  cloud  standard  deviation  in  the  x  and  y  directions. 
Physically,  Hopkins*  smearing  equation  calculates  the  dose  from  a 


Gaussian  cloud  for  points  which  are  on  a  line  perpendicular  to  the 
effective  wind  vector.  Or,  given  the  effective  wind  to  a  point  on  the 
hotline,  Eq(l)  calculates  the  dose  rate  for  points  in  the  crosswind 
direction.  Appendix  B  has  an  explanation  of  the  geometry  of  the  smearing 
equation  and  its  physical  significance. 

The  source  normalization  constant  typically  has  values  of  2000  to 
2950  R-mi  sq/hr-kt(4 :212; 3 : 1 ) .  The  dose  program  uses  2950  R-mi  sq/hr-kt 
or  7452  rads-tissue-km  sq/hr-kt. 

The  normalized  fractional  rate  of  arrival  of  radioactivity  at  the 


time  of  arrival  is  calculated  from  a  method  developed  by  Bridgman  and 
Bigelow(4:210).  The  rate  of  arrival  is: 


-  the  activity  size  distribution,  or  the  fraction  of  activity 
in  particles  of  size  r  (per  micron) 

*  the  rate  of  change  of  particle  radius  arriving  on  the 
g round ( mi c rons / hour ) 

The  activity  size  distribution  is  assumed  to  be  log  normal,  and  r  and 
dr/dt  are  calculated  from  a  Laurent  series  fit  using  coefficients 
calculated  by  Colarco(4:214;7:67). 


Dose  at  Arbitrary  Coordinates 

The  hotline  position  is  known  from  the  hotline  locator's  output. 

But  to  find  the  dose  at  a  cell  center,  it  is  necessary  to  find  the 
hotline  coordinates  such  that  the  cell  center  is  in  the  crosswind 
direction.  Figure  3  shows  a  hypothetical  hotline,  made  up  of  hotline 
segments,  and  a  cell  center.  The  equation  of  the  dotted  line  that 
contains  the  second  hotline  segment  can  be  determined  from  the  trace 
particle  coordinates.  It  turns  out  that  there  are  two  points  on  the 
dotted  line  such  that  the  cell  center  is  in  the  crosswind  direction.  If 
either  of  these  points  is  on  the  hotline  segment,  then  those  hotline 
coordinates  are  used  in  Eq(l)  to  calculate  the  dose  at  cell  center.  If 
both  points  are  on  the  hotline  segment,  then  the  one  with  the  earliest 
time  of  arrival  is  used  to  calculate  the  dose. 

Appendix  C  has  a  derivation  of  the  quadratic  equation  used  to  find 
the  hotline  coordinates,  given  a  cell  center  location.  The  derivation 
also  explains  how  the  time  of  arrival  and  wind  shears  are  calculated  for 


hotline  points  between  the  trace  particle  positions. 
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FIGURE  3  Hotline  Position  (X,Y)  and  Cell  Center  (x,y)  Geometry 


Population  Model  and  Mesh  Size 

The  population  model  was  developed  by  Wendel(21:5-7).  In  the 
Continental  US  there  are  316  Standard  Metropolitan  Statistical 
Areas(SMSA),  which  contain  the  urban  population  of  the  United  States. 
Wendel  found  the  approximate  geographic  centroid  of  each  SMSA  and  placed 
the  SMSA  population  at  the  centroid.  Approximately  one  quarter  of  the 
population  lives  in  rural  areas.  Using  the  rural  population  density  of 
each  state,  Wendel  was  able  to  calculate  a  homogeneous  (state  by  state) 
rural  population  within  each  one  degree  latitude  by  one  degree  longitude 
cell.  Wendel  then  assumed  that  the  population  of  each  cell  resided  at 
its  center. 


AFOTEC  wanted  to  know  the  dose  rate  for  a  mesh  size  of  about  50 
miles  by  50  miles  across  the  US.  A  one  degree  by  one  degree  mesh  is 
about  69  miles  by  55  miles,  so  this  mesh  size  satisfied  AFOTEC 
requirements  and  was  a  convenient  mesh  size  for  calculating  population 
dose.  However,  dose  rates  for  this  mesh  size  are  sensitive  to  minor 
changes  in  burst  position.  That  is,  if  the  dose  rate  at  a  given  cell 
center  kills  zero  percent  of  the  population,  then  moving  the  burst  point 
ten  or  15  kilometers  could  give  a  dose  rate  high  enough  to  kill  100 
percent  of  the  cell  population.  Similarly,  by  assuming  that  the  entire 
population  of  an  SMSA  resides  at  its  centroid,  a  small  change  in  burst 
position  has  a  large  effect  on  that  population's  dose. 

To  correct  this  problem,  the  dose  program  calculates  the  average 
dose  across  a  cell  or  SMSA,  and  assumes  that  the  population  is  evenly 
distributed  across  each  cell  or  SMSA.  The  average  dose  can  be  found  from 
a  difference  of  cumulative  normal  functions.  (See  Appendix  D)  The 
average  SMSA  has  about  4600  square  kilometers(19:26) ,  and  the  dose 
program  assumes  that  each  SMSA  is  a  circle  with  a  radius  of  38 
kilometers . 

Although  using  the  average  dose  allows  the  use  of  a  coarse  mesh  for 
the  population,  the  results  can  be  deceiving.  Consider  Figure  4,  which 
shows  a  side  view  of  the  dose  deposited  across  a  cell,  on  a  line  through 
the  cell  center.  Let  the  fallout  be  from  a  one  megaton  burst  with  500 
kilotons  fission;  let  the  effective  wind  be  50  kilometers  per  hour;  let 
the  time  of  arrival  be  three  hours  and  the  cloud  standard  deviation  in 
the  crosswind  direction  be  six  kilometers.  Under  the  hotline,  the  time 
integrated  dose  to  infinity  is  1030  rads  tissue.  For  someone  living  more 
than  three  standard  deviations  from  the  hotline,  the  dose  is  zero.  The 
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FIGURE  4  Dose  Across  a  Cell 

average  dose  across  the  cell  is  140  rads  tissue.  In  this  case,  the 
average  dose  is  almost  eight  times  less  than  some  people  receive.  And 
the  average  dose  predicts  no  casualties  while  over  98  percent  of  the 
people  directly  below  the  hotline  will  actually  die.  If  several  hotlines 
cross  the  cell,  then  the  average  dose  will  give  a  more  accurate  account 
of  the  deaths  within  the  cell.  Also,  at  later  times  when  the  cloud  size 
is  on  the  order  of  or  larger  than  the  cell  size,  the  average  dose  will 
more  accurately  predict  population  deaths.  There  was  insufficient  time 
to  do  a  complete  analysis  of  the  effects  of  using  the  average  dose. 

Multiple  Bursts 

When  there  is  an  attack  on  an  area  target,  such  as  a  Minuteman 
missile  field,  it  is  desirable  to  treat  the  fallout  produced  as  one 
mega-cloud.  This  allows  calculating  the  dose  rate  for  a  cell  one  time. 


V  V  *  .-  .*  *„•  V  " 
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instead  of  superimposing  the  dose  rates  from  several  hundred  bursts. 
Crandley  developed  and  validated  a  multiple  burst  technique  for  a 
constant  wind  fallout  code(8).  Appendix  E  draws  on  Crandley' s  results 
and  extends  them  to  variable  winds.  There  was  insufficient  time  to 
derive  and  validate  expressions  for  the  average  dose  across  a  cell  or 
city  from  an  area  of  bursts.  Therefore,  as  will  be  shown  later,  the  area 
of  bursts  must  be  on  the  order  of  or  larger  than  the  cell  size  to  achieve 
consistent  dose  rate  calculations. 

Probability  of  Death 

The  probability  of  death  as  a  function  of  gamma  radiation  dose 
received,  can  be  inferred  to  be  a  cumulative  log  normal  distribution  (See 
Figure  5)  with:  170  rads  tissue  causing  a  one  percent  chance  of  death, 
1200  rads  tissue  causing  a  99  percent  chance  of  death,  and  450  rads 
tissue  causing  a  50  percent  chance  of  death  within  60  days( 15 : 20-22) . 

This  probability  distribution  is  valid  for  a  dose  received  over  a  short 
time  period  of  a  few  hours.  The  Way-Wigner  decay  formula  is  used  to  find 
the  time  integrated  dose  from  the  fallout  time  of  arrival  to 
inf inity(20: 1318).  It  was  assumed  that  the  dose  to  infinity  was 
delivered  in  a  short  time  period.  The  effects  of  alpha  and  beta 
radiation,  and  the  effects  of  fallout  ingestion  were  ignored. 

Glasstone  gives  typical  values  for  the  dose  protection  factor 
afforded  by  various  types  of  shelters(10:441) .  A  protection  factor  of 
three  was  assumed  to  be  a  nominal  value  for  the  rural  and  urban 
population.  The  dose  program  uses  the  log  normal  distribution  above,  and 
a  dose  protection  factor  of  three  to  calculate  the  probability  of  death. 
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Wind  Shear  and  Cloud  Size 

Wind  shear  acting  on  the  fallout  cloud  as  it  translates  through  the 
atmosphere  will  cause  the  cloud  to  spread(12 : 49 ) .  The  AFIT  Smear  Model 
uses  the  WSEG-10  formulation  for  the  cloud  standard  deviation  as  a 
function  of  time(3:2): 

r-  fr*  *  (GS*)*]*  to) 

where 

*  the  cloud  standard  deviation  in  the  x  or  y  direction 
A?  =  the  time  of  arrival 

*  the  cloud's  initial  standard  deviation  in  the  vertical 
S  ~  the  wind  shear  in  the  x  or  y  direction 

Or  -  a  term  for  toroidal  growth  of  the  cloud  at  early  times. 

After  three  hours,  for  a  one  megaton  burst,  Qy  is  a 
constant  six  kilometers. 

Typical  values  of  wind  shear  used  are  on  the  order  of  one  to  two  per 
hour,  which  gives  a  cloud  standard  deviation  of  50  to  100  kilometers  at 
late  times  for  a  one  megaton  burst.  The  hotline  locator  program 
calculates  wind  shear  values  of  about  one  to  30  per  hour,  which  in 
general  agrees  with  the  range  of  wind  shear  calculated  by  Nerment  for 
certain  selected  wind  fields(17: 57) .  Typical  values  of  the  hotline 
locator  wind  shear  at  late  times  for  a  one  megaton  burst,  give  cloud 
standard  deviations  of  about  500  to  1000  kilometers.  But  Bauer  has 
determined  that  the  cloud  standard  deviation  can  be  calculated 


from(18:Sec  6,39): 
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where 

=  a  latitude,  altitude  and  time  of  year  dependent  variable 
7*  «  the  time  after  burst. 

With  QT  in  kilometers,  T  in  days  and  using  a  mean  value  for 

res  300  rty* 


(6>) 


So  at  late  times,  the  cloud  standard  deviation  should  be  on  the  order  of 
200  kilometers. 

Norment  has  developed  a  smear  fallout  code  called  DNAF-1.  His 
empirical  equation  for  horizontal  cloud  standard  deviation  is(17:41) 

<r=  [%* y  (AiS-tzA a)*] *  / 7) 


with 

£ 


o  ■  a  term  for  growth  not  due  to  wind  shear 
»  the  vertical  thickness  of  the  cloud 

ta 

m  the  time  of  arrival 

3  *  the  wind  shear. 


The  dose  program,  following  Norment,  uses: 

r - 


WgfcMU*  |  HV1.IUV11V  | 
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to  calculate  the  standard  deviation.  This  "empirical  fix”  of  the  WSEG-10 
formula  gives  cloud  standard  deviations  of  about  100  to  200  kilometers 
for  typical  values  of  wind  shear  and  a  24  hour  time  of  arrival,  which 
approximates  Bauer's  prediction. 
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Program  Validation 

The  dose  program  was  used  with  a  variety  of  hypothetical  hotlines  to 
check  its  operation.  For  instance,  reflecting  the  hotline  about  the  x  or 
y  axis  results  in  the  dose  rates  being  reflected  about  the  x  or  y  axis. 

As  a  further  test,  a  conservation  check  was  made  of  the  activity  smeared 
on  the  ground.  For  single  and  multiple  burst  hotlines,  the  burst 
position  was  varied  within  a  latitude  longitude  cell,  and  the  program 
output  was  summed  to  find  the  total  activity  on  the  ground,  as  compared 
to  the  activity  contained  in  20  micron  and  larger  particles. 

When  activity  conservation  was  checked  for  a  single  burst  with  a  one 
degree  by  one  degree  cell  size,  and  dose  averaging  was  not  used,  then 
anywhere  from  20  to  180  percent  of  the  available  activity  was  on  the 
ground.  With  dose  averaging,  101  +  13  percent (3  sigma  deviation,  20 
locations)  of  the  activity  was  smeared  on  the  ground.  Without  dose 
averaging,  a  cell  size  of  ten  km  by  ten  km  conserved  activity,  and  the 
run  time  for  one  hotline  increased  from  five  to  50  seconds. 

For  multiple  burst  hotlines,  dose  conservation  varied  depending  on 
the  size  of  the  burst  area.  Table  II  shows  how  dose  conservation  varied. 
It  is  recommended  that  an  area  of  bursts  have  dimensions  of  at  least  100 
by  100  kilometers.  Otherwise,  the  bursts  should  be  treated  by 
superposition. 

Results 


For  an  attack  against  the  25  air  bases  in  Table  I,  the  dose  program 
used  30  seconds  of  central  processor  time  and  calculated  537,000  fallout 
deaths  within  60  days  of  the  attack.  No  single  burst,  by  itself,  had  a 


TABLE  II 


Activity  Conservation  for  Multiple  Bursts 
Versus  Size  of  Burst  Area 


Number 

Area 

X  Activity  on 

of  Bursts 

Dimensions(km) 

the  Ground 

100 

50 

X 

50 

233(+  121) 

100 

100 

X 

100 

135(  +  38) 

165 

156 

X 

180 

112(  +  10) 

220 

180 

X 

208 

108(  +  11) 

100 

200 

X 

200 

102(  +  13) 

Deviations  are  three  sigma  for  ten  locations.  Point  cell 
values  are  used.  That  is,  dose  is  not  averaged  across  a  cell. 


dose  rate  high  enough  to  cause  any  deaths,  due  to  the  effects  of 
averaging  the  dose  across  a  cell.  So  the  deaths  occurred  in  the  areas 
where  several  hotlines  overlap.  (See  Figure  2)  For  an  attack  using  a  one 
megaton  surface  burst  against  each  Minuteman  silo  and  control  center,  the 
dose  program  used  13  seconds  and  calculated  25.9  million  deaths.  A 


combined  attack  used  37  seconds  and  calculated  26.9  million  deaths. 


The  hotline  locator  and  dose  programs  use  a  reasonable  amount  of 
computer  time  to  calculate  the  dose  and  population  insult  across  the 
United  States.  To  calculate  the  hotlines  and  the  dose  from  the 
hypothetical  attack  with  31  hotlines  required  407  central  processor 
seconds,  on  a  CYBER  845.  The  technique  of  averaging  the  dose  across  a 
grid  cell,  allows  the  use  of  a  large  mesh  and  results  in  a  faster  run 
time  but  underestimates  the  fallout  deaths  for  single  bursts.  Several 
hundred  bursts  in  one  area  can  be  treated  as  a  mega-cloud  with  one 
hotline.  However,  the  dimensions  of  the  area  should  be  as  large  or 
larger  than  the  cell  size  to  conserve  activity  on  the  ground. 

There  are  several  areas  in  which  the  models  or  programs  can  be 
improved : 

1.  An  investigation  of  the  growth  of  fallout  clouds  due  to  wind 
shear  needs  to  be  performed  to  derive  an  accurate  expression  for 
cloud  standard  deviation. 

2.  The  spectral  coefficients  used  to  calculate  wind  speed  can  also 
be  used  to  calculate  vertical  wind  speed  in  the  atmosphere(13) .  The 
vertical  winds  should  be  added  to  the  hotline  locator  program. 

Since  Colarco's  Laurent  coefficients  for  r(t)  were  developed 
assuming  no  vertical  winds,  this  would  require  a  different  method  be 
used  to  calculate  g(t). 

3.  An  algorithm  to  find  the  average  dose  across  a  cell  for  multiple 
bursts  should  be  added  to  the  dose  program. 


4.  When  the  fallout  cloud  dimensions  are  larger  than  the  synoptic 
scale  of  the  amosphere  (about  150  to  450  km),  a  cloud  may  tend  to 
disperse  into  several  smaller  clouds  which  translate  separately( 13) . 
An  investigation  should  be  done  to  determine  how  local  fallout  is 
affected  by  synoptic  scale  wind  patterns. 

5.  The  hotline  locator  program  should  be  further  modified  to  reduce 
its  run  time.  The  spectral  wind  coefficients  for  different  times 
could  then  be  interpolated  to  find  hotlines  in  a  spatially  and 
temporally  varying  wind  field,  with  a  reasonable  run  time. 


Figures  6-11  show  the  fallout  cloud  locations  from  the  counterforce 
attack  in  Table  I.  The  marker  locations  show  where  the  fallout  clouds 
are  touching  the  ground  for  every  four  hours  after  the  attack.  The 
markers  show  only  where  the  clouds  are  touching  the  ground.  Wind  shear 
will  cause  the  cloud  position  to  change  with  altitude.  The  solid  lines 
show  where  the  fallout  clouds  have  already  passed,  and  the  dotted  lines 
show  thw  future  track  of  the  clouds. 
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Appendix  B 

Physical  Interpretation  of  Hopkins  *  Variable 
Wind  Smearing  Equation 


The  physical  significance  of  Hopkins'  variable  wind  smearing 
equation  was  first  shown  by  Bridgman(2).  This  discussion  generally 
follows  his  explanation. 

Consider  Hopkins'  variable  wind  smearing  equation  written  in  terms 
of  velocities,  rather  than  hotline  coordlnates(ll :83) : 


where 

D>H)  -  the  Unit  Time  Reference  Dose  Rate  at  (x,y) 

K  ■  the  source  normalization  constant 

Xf  -  the  fission  yield 

■  the  fractional  rate  of  arrival  of  radioactivity  at 
the  time  of  arrival 

Vifc  ■  the  effective  wind  velocity  in  the  x  direction,  or  x/t& 

the  effective  wind  velocity  in  the  y  direction,  or  y/te, 

the  cloud  standard  deviations  in  the  x,y  directions 
Figure  12  shows  the  elliptical  two  dimensional  Gaussian  cloud 
located  at  (X,Y).  How,  define  the  effective  wind  vector  as: 


v&  =  vfc  £  +  Vff 

therefore 

Ve*  I Vel'(V*+  Vff* 


(*’*) 

(3-3) 
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Figure  12 

Two  Dimensional  Gaussian  Fallout  Cloud 


Define  the  effective  wind  vector  angle  as: 

<?(  =  /(y/x)  (tf-y) 

Then,  the  cloud  standard  deviation  vector  in  the  crosswind  direction  is 


r 


—  i* 


so  that 
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Now,  from  Eq(B-l)  consider  the  term 
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Now,  consider  the  exponential  term  in  Eq(B-l) 
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From  Figure  12,  the  unit  vector**,  in  the  crosswind  direction  from  (x,y) 
to  (X,Y)  is: 


jy  -  £ 


/ a~  ui ) 


effective  wind. 


with  the  wind  shear  and  time  of  arrival  for  each  trace  particle,  this 


method  will  find:  the  hotline  coordinates  which  define  the  effective  wind 
vector,  the  time  of  arrival  and  the  wind  shear. 

Figure  13  shows  the  ith  hotline  segment  of  a  hotline,  a  known  point 


(x,y),  and  the  unknown  hotline  coordinates  (X,Y). 


0 - o 

HOTLINE  SEGMENT 


The  equation  of  the  line  which  contains  the  ith  hotline  segment  is 

/  =  afr)X  (C-t) 


where 
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and 

b(i)  -  yu)  -  Xh)  (c-3 

Defining  the  vectors  R,  P,  r  as: 

A=X^"  y/  (c-4) 

?■=■*£+?$  (c-f) 

p-(- *-x)£  y)/  (c-t,) 


Then,  since  (x,y)  is  in  the  crosswind  direction, 

4-P  -  O 


(c-r) 


(y-y) x  +  (?-y) y  (c-s) 


Substituting  Eq(C-l)  for  Y  into  Eq(C-8)  gives 


AX1  +  6X.  +  C*# 
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Therefore , 
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Eq(C-ll)  gives  two  trial  Xs  on  the  straight  line  containing  the  ith 
hotline  segment,  such  that  (x,y)  is  in  the  crosswind  direction.  If 
neither  trial  X  is  on  the  ith  segment,  then  the  ith  plus  one  segment  must 
be  checked.  If  only  one  trial  x  is  on  the  ith  segment,  then  that  (X,Y) 
defines  the  hotline  coordinate.  If  both  trial  Xs  are  on  the  hotline 
segment,  then  the  one  with  the  earliest  time  of  arrival  defines  (X,Y). 

To  find  the  time  of  arrival  and  wind  shear  associated  with  the 


hotline  coordinate  (X,Y),  let's  assume  that  the  time  of  arrival  and  wind 
shear  vary  linearly  along  the  ith  segment.  Then  with: 


ttL=  T±a.(iri)  + [l-rjta.fr)  ( 

S*  -  TUti*/)  +U-T]S*fr)  (C 

5*  =  TS^^t)  +D~T]Sffr)  (t 


with 


t<L 


the  time  of  arrival 


the  wind  shear  in  the  x,y  directions 


Now  these  hotline  variables  can  be  used  to  calculate  the  unit  time 


reference  dose  rate  at  (x,y). 


Appendix  D 


Average  Dose  Across  £  Cell 


In  Figure  14,  the  dose  rate  at  (x,y)  is: 

¥  y^x)3- 


D.faf)  ~ 


\[#r  r  £ 


(o-O 


with 


C-  k Yf  f(ta) to. 

<r-  (t*y*+fs.*xx)* 

and  with  (X,Y)  on  the  hotline.  The  average  dose  rate  along  the  line 


is(14:739) : 


4**  = 


(D-A) 


At  later  times  when  fty  changes  little  in  the  time  it  takes  for  the 
fallout  cloud  to  cross  the  cell,  then  the  average  dose  rate  calculated  by 
Eq(D-2)  is  the  average  dose  rate  across  the  entire  cell. 

Defining  as: 


c*  ■=  Y 

(o-3) 

line  of  length 

'jL  ^  X  -'f' ^ 

(p-4a) 

(0-46) 

when  $ 


is  zero  on  the  hotline. 


Figure  14 

Average  Dose  Geometry 


Substituting  Eqs(D-4)  into  Eq(D-2)  gives: 


•  __c_ 

b -a 


\j=Tr 


2  "  r 


(~y44*''°(~'XC04'Ot ) 


and  substituting  (D-6)  and  (D-7)  into  Eq(D-5)  gives  the  average  dose  rate 


as : 

4 

0,44 


£ 

(6  -t 2 


where 


T 


Cos) 

(0-4) 


and 


~S  Y4+*q(  -  4 

r 


The  integral  in  Eq(D-8)  is  the  cumulative  normal  function,  which  can  be 
calculated  by  a  polynomial  approximation  f.:om  Abramowitz  and 
Stegun(l :932) . 

To  find  the  line  lengths  a  and  b,  definefl^^y  as: 


^^2/  =  '&*?*(&/#) 


where 


//  =  the  north-south  dimension  of  the  cell 
Id  53  the  east-west  dimension  of  the  cell. 


Then  with  d  defined  as: 


j,L(*.-x)\(ry)] 


</* 


Co-//) 


(o-u) 


when OismAjC  then 


Appendix  E 


Dose  From  Multiple  Bursts  Within  a  Rectangular  Area 

The  following  method  finds  the  dose  rate  at  a  cell  center  from  an 
area  of  bursts.  Three  primary  assumptions  are  made.  First,  the  fallout 
observer  must  be  far  downwind  compared  to  the  area  dimensions,  so  that 
the  fallout  time  of  arrival  does  not  vary  significantly  in  the  amount  of 
time  it  takes  the  mega-cloud  to  pass  his  position.  This  means  that 

does  not  change  significantly  for  each  burst,  and  the  observer 
cannot  tell  whether  the  fallout  was  deposited  by  an  area  or  a  line  of 
bursts.  The  second  assumption  is  that  the  summation  of  dose  rates  over 
the  line  of  bursts  can  be  replaced  by  a  line  integral.  Crandley  showed 
that  this  assumption  is  valid  when  the  burst  separation  is  less  than  one 
cloud  standard  deviation(8:5-ll) .  Finally,  it  is  assumed  that  each  cloud 
in  the  area  of  bursts  is  translated  by  the  wind  affecting  the  center 
burst.  The  development  below  finds  the  burst  density  along  a  crosswind 
line  through  the  center  of  the  area  of  bursts  then  uses  the  burst  density 
to  calculate  the  dose. 

Variable  Burst  Density  Geometry 

Consider  the  area  of  bursts  in  Figure  15,  with  height  H,  width  W  and 
with  the  bursts  evenly  distributed  over  the  area.  An  observer  far 
downwind  cannot  tell  whether  the  fallout  was  deposited  by  the  area  of 
bursts  or  by  a  line  of  bursts  in  the  crosswind  direction. 

Define  as  anS^e  between  the  x  axis  and  the  rectangle's 


diagonal,  or 


Figure  16  shows  that  when  o(  =  0,  or  the  effective  wind  is  in  the 


x  direction,  that  the  burst  density  is  a  constant  in  the  crosswind 
direction.  In  the  more  general  case  when  O  4  of <  •  Figure  17  shows 

that  the  burst  density  is  a  trapezoid.  The  area  of  the  trapezoid  is: 

h  (£- 3 j 

which  is  also^'th^total  number  of  bursts  within  the  target  area.  The 
diagonal  of  the  targe t "irre^i s  : 

O/A  -  (£-4) 

and  from  figure  17,  the  base  of  the  trapezoid,  which  is  the  crosswind 
line  of  length  $  is: 

./=  Of  A  to+P  (£-$) 

where  j3  is  the  angle  from  the  crosswind  line  to  the  diagonal  in  right 
triangle  ABC. 

Notice  that 

so  that 

J.  £/A  (£-7) 

Now,  define  AP  as  the  line  length  DE,  then  the  distance  EF  is 


equal  to  AG,  or 


BURSTS/KM 


*. 

I  L 


Figure  18 

Trapezoidal  Burst  Density  Geometry 
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From  right  triangle  ADE, 


AP  = 


(£- 


and  from  right  triangle  AGH  and  Eq(E-8), 


b/£_ 

H/H-AP 


(£-/ 


where  b  is  the  length  of  the  top  of  the  trapezoid  in  Figure  17. 
Substituting  Eq(E-9)  into  Eq(E-lO)  and  solving  for  b  yields: 


b  -  (M M,* 


Since  the  area  of  the  trapezoid  is  equal  to  the  total  number  of 


bursts  N,  then  from  Eq(E-3), 


AAl_ 


where  b  is  known  from  Eq(E-ll)  and  is  known  from  Eq(E-7). 

In  a  similar  manner  to  the  above ,  whenp^.^^  then: 


J  -  A/P  ( t****.)  (£ 


Dose  From  a  Line  of  Bursts 


Crandley  has  shown  that  when 


c  a/  r  / 

D'(*'t)=  T  * 

Lth 

C  -  A 

X  =  the  x  value  of  the  hotline  coordinate 


0,  as  in  Figure  16,  that (8: 8), 

L 


~r‘dS  (£-&) 


vy^=  the  cloud  standard  deviation  in  the  crosswind  or  y  direction 

Af  =  the  number  of  bursts 

//  =  the  length  of  the  line  of  bursts 


H/H  - 


the  constant  burst  density  along  the  line  of  bursts 


Figure  19  shows  a  line  of  bursts  of  length 


in  the  crosswind 


direction  depositing  activity  at 


(**/)■ 


Extending  Crandley' s  results 


to  a  variable  burst  density  and  to  the  variable  wind  smearing  equation. 


Eq(E-16)  becomes 


U/i 


tit  . 

gag) 

.a fafrr  e 


JS 


where 


the  burst  density  as  a  function  of  the  distance  S  along 


the  line  of  bursts 

/ \/i  .  /t-.  *  ^ 


r  ■( 

and  vt  are  functions  of  S. 


By  definition,  each  burst  point  is  the  coordinate  origin  for  the 


fallout  from  that  burst(ll :73-82) .  Since  the  same  wind  is  assumed  to 


Figure  19 

Dose  From  a  Line  of  Bursts 


translate  all  clouds,  notice  in  Figure  19  that  every  burst  on  the  line 


has  the  same  effective  wind  vector  to  its  hotline,  and  that  X  and  Y  are 


the  same  for  each  burst, 
different  for  every  burst 


However,  the  crosswind  coordinates 
Redefining  as: 


and 

tri  -  ^ 

(yj^+XCt**) 

For  BD(S)  in  Eq(E-20),  consider  the  trapezoidal  burst  density  in 
Figure  16. 


and  when  ^  ±  Tp 


3Q(s)  - 


-IAS 

w 


(£-  32  c) 


Substituting  Eqs(E-21)  and  Eqs(E-22)  into  Eq(E-20),  we  see  that  the 
integral  in  Eq(E-20)  can  be  broken  up  into  three  integrals.  The  first 


integral  is: 


3h  f  n-fry-fXi  Jjh  7 _ I  "sM  j  f£-13a,} 


UL  - 


LL- 


r 

'±  -  £-U~«)y-(jl.+  ^rCatet)X  (£-246, 


The  second  integral  is: 


JpUL* 

lu  i 


UL> ,  f*j4+^)Y-/f-4*v<)x_ 

r 

/,<  fr-4-om*t)y-fe+4-4>**<2X 


(£-336) 

(E-2fd.) 


/iT-  ) 


where  CNF  is  the  cumulative  normal  function,  which  can  be  calculated  by  a 
polynomial  approximation  from  Abramowitz  and  Stegun(l :932) . 

Eq(E-23b)  reduces  to 

=  —  [CUF(UL')  -fA/FCU')]  (S-JU) 

Comparing  Eq(E-23a)  with  Eq(E-23c),  and  comparing  Eqs(E-24)  with 
Eqs(E-26),  we  can  write  down  the  results  of  Eq(E-23c)  directly  as: 


Now,  the  dose  rate  from  an  area  of  bursts  can  be  found  from: 

4  (<f)  =  (£'au)  +  bte-in)  (E-tf) 


Appendix  F 

Hotline  Locator  Program  Listing 

This  program  Is  the  modified  version  of  Hopkins'  hotline  locator 
described  In  Section  II. 

Table  III  is  an  example  Tape20,  a  burst  data  file.  The  first  line 
is  for  a  single  burst,  and  the  second  line  is  for  multiple  bursts  in  an 
area  target.  The  entries  in  order  are:  center  burst  latitude,  center 
burst  longitude,  weapon  yield  in  megatons,  number  of  bursts,  height  of 
target  and  width  of  target  in  kilometers,  the  yield  fission  fraction,  the 
activity  size  distribution  volume  fractionation,  natural  log  of  the 
activity  size  distribution's  mean  radius t  activity  size  distribution 
logarithmic  slope,  and  particle  density  in  grams  per  cubic  centimeter. 

The  activity  size  distribution  is  assumed  to  be  log  normal  as  described 
by  Bridgman  and  Bigelow(4:210-211).  The  values  in  Table  III  are  for  the 
DELFEC  default  distribution.  The  hotline  locator  will  work  with  any 
particle  density,  however,  the  dose  program  will  only  work  with  a 
particle  density  of  2.6  grams  per  cubic  centimeter. 

Table  IV  is  an  example  Tape21,  the  program's  output  file.  The  first 
line  is  the  burst  data  for  the  hotline.  Following  that  is  one  line  of 
data  for  each  trace  particle,  showing  in  order:  trace  particle  radius  in 
meters,  time  of  fall  in  seconds,  x  position  in  meters,  x  wind  shear  in 
per  second,  y  position  in  meters,  and  y  wind  shear  in  per  second. 
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TABLE  III 


Example  Tape20  Burst  Data  File 


47.95  97.40 

47.50  110.00 


.  1Q0E+01  1 

. 100E+01  220 


000  000  .50 

180  208  .50 


.68  -1.58964 

.68  -1.58964 


1.38629  2.6 

1.38629  2.6 


TABLE  IV 


Example  Tape21  Hotline  Data  File 


47.95  97.40 

. 100E+01 

1.  0.  0. 

.50  .68 

-1.58964 

1.38629 

. 200E-04 

. 942E+05 

. 176E+07 

. 316E-02 

.  440E+06 

. 393E-02 

.213E-04 

.840E+05 

. 169E+07 

.310E-02 

. 502E+06 

. 389E-02 

.225E-04 

.754E+05 

. 163E+07 

.304E-02 

.539E+06 

.385E-02 

. 250E-04 

. 621E+05 

.  151E+07 

. 291E-02 

. 552E+06 

. 374E-02 

.275E-04 

. 522E+05 

. 136E+07 

. 282E-02 

.511E+06 

. 363E-02 

. 300E-04 

. 448E+05 

. 121E+07 

•285E-02 

. 441E+06 

. 358E-02 

. 338E-04 

.365E+05 

. 102E+07 

. 312E-02 

.326E+06 

.  34  9E-02 

. 400E-04 

.275E+05 

.779E+06 

.366E-02 

. 181E+06 

.  337E-02 

. 500E-04 

. 178E+05 

. 505E+06 

. 432E-02 

. 527E+05 

. 331E-02 

. 875E-04 

. 821E+04 

.238E+06 

.444E-02 

- . 704E+04 

.  267E-02 

47.50  110.00 

. 100E+01 

220.  180.  208. 

.50  .68 

-1.58964 

1.38629  : 

. 200E-04 

. 942E+05 

.265E+07 

. 287E-02 

- . 788E+06 

.581E-02 

. 213E-04 

. 840E+05 

.232E+07 

. 273E-02 

-.797E+06 

.547E-02 

. 225E-04 

. 754E+05 

.  206E-1-07 

•272E-02 

-.905E+06 

. 520E-02 

. 250E-04 

.621E+05 

.  163E+07 

. 309E-02 

- . 108E+07 

.485E-02 

. 275E-04 

. 522E+05 

. 137E+07 

.351E-02 

-. 111E+07 

.473E-02 

. 300E-04 

. 448E+05 

. 120E+07 

•380E-02 

-.105E+07 

.473E-02 

. 338E-04 

.365E+05 

. 104E+07 

.408E-02 

-. 922E+06 

.483E-02 

. 400E-04 

. 275E+05 

. 844E+06 

. 433E-02 

- . 747E+06 

.492E-02 

. 500E-04 

.  178E+05 

.565E+06 

.422E-02 

-.522E+06 

.488E-02 

.875E-04 

. 821E+04 

. 251E+06 

.  386E-02 

- . 259E+06 

.475E-02 

PROGRAM  HBATCH 

C  THIS  IS  A  TRIAL  RUN  OF  HBATCHN  CALLED  H0LINE7  23  NOV  84 
C  MODIFIED  FOR  TAPE20  INPUT 

C  THIS  PROGRAM  COMPUTES  THE  TRANSLATION  OF  PARTICLES 
C  FALLING  FROM  THE  INITIAL  CLOUD  THROUGH  SPECTRAL  WINDS. 

C  THE  ATMOSPHERE  IS  DISCRETIZED  INTO  NL  LAYERS,  FROM  EARTH 
C  SURFACE  UP  TO  THE  HIGHEST  PARTICLE. 

C  WIND  VECTORS  ARE  LINEARLY  INTERPOLATED  INTO  A  LAYER  BY 
C  INTERPOLATING  WINDS  FROM  SPECTRAL  LEVELS  ABOVE  AND  BELOW. 

PARAMETER  (NL=12,NP=10, JCAP-30, PI=3. 14159265) 

C 

C  NL*  NUMBER  OF  LAYERS  IN  ATMOSPHERE 
C  NP  =  NUMBER  OF  TRACE  PARTICLES  USED  TO  DEFINE  HOTLINE 
C  CHOOSE  DELZ  GT  100  METERS  (BL  THICKNESS) 

C 

DIMENSION  PDIA(NP) , ZPDIA(NP) ,WT(NP) 

DIMENSION  ZMSL (NL+1 ) , HGEO (NL+l ) , TZ (NL+1 ) , PR (NL+1 ) , 

C  DENS (NL+1 ),VISK (NL+1) 

DIMENSION  TR (NL) 

DOUBLE  PRECISION  CC, COL RAD 

DOUBLE  PRECISION  RR,RADLON, EPS (JCAP+2, JCAP+1) 

COMPLEX  CU ( 1 s 12, 1 ! JCAP+2, 1 s  JCAP+1 ) , CV ( 1 s 12, 1 : JCAP+2, 1 : JCAP+1 ) 

REAL  BLAT, BLON, YM, NBURS, HLINE, WLINE, FF, FV, PALPH, PBETA, PDEN 
C 

CALL  EPSLN (EPS, JCAP) 

REWIND  10 
REWIND  20 
REWIND  21 
DO  99  LVL-1,12 

READ ( 10, 98) ( (CU (LVL,  N, L) , N*1 , 32) , L-l , 31 ) 

READ (10, 98) ( (CV(LVL,N, L) , N«1 , 32) , L=1 , 31 ) 

99  CONTINUE 
98  FORMAT (6E13. 7, 2X) 

C 

C 

29  READ (20,30, END-3 1 >  BL AT , BLON , YM , NBURS , HL I NE , WL I NE , FF , FV, 

C  PALPH, PBETA, PDEN 

WR I TE ( 2 1 , 30 ) BL AT , BLON , YM , NBURS , HL I NE , WL I NE , FF , F V , 

C  PALPH, PBETA, PDEN 

30  FORMAT (F5. 2, F8. 2, El 1.3, 3F5. 0, F5. 2, F6. 2, 2F10.5,F5. 1) 

COLRAD-P 1/2. -BLATtPI / 1 80 . 

RADLON-2 .  tP  I -BLON  *P I / 1 80 . 

CC-COLRAD 

RR-RADLON 

C 

YKT-YM41000. 

DATA  PDIA/.004, . 00425, . 0045, . 005, . 0055, . 006, . 00675, . 008, .  01, .0175/ 
C 

C  COMPUTE  TRACE  PARTICLE  HEIGHTS  IN  THE  STABILIZED  CLOUD 


onnononnn  n  no 


USING  DELFIC  CORRELATIONS  IN  H0S2 


YL-LOG(YKT) 

SL-1.574-.01197*YL+.03636*YL*YL-.0041*YL**3+.0001965*YL**4 

SLOPE«EXP(SL> 

X I -7 . 889+ . 34* YL+ . 00 1 226* YL* YL- . 005227 * YL*  *3+ . 0004 1 7*YL*  *4 
XINT-EXP(XI) 

DO  1  1*1, NP 

ZPDI A ( I ) *-SLOPE*PDIA ( I ) * 10000.  +  XINT 
1  CONTINUE 

ZMAX  IS  THE  HEIGHT  OF  THE  SMALLEST  TRACE  PARTICLE: 

THE  MAX  HEIGHT  OF  WAFER  CENTERS  IN  THE  CLOUD 
ZMAX*ZPDIA(1) 

DIVIDE  THE  ATMOSPHERE  INTO  NL  LAYERS  FROM  0  TO  ZMAX 
DELZ*ZMAX/NL 

COMPUTE  ATMOSPHERE  STATE  PROPERTIES 

CONSTANTS  IN  ATMOSPHERE  STATE  EQUATIONS 
RERTH*6356766 . 

GAMMA- 1. 

GOP-9. 80665 

XM0-28.9644 

RSTAR-8314.32 

BETA-1.458E-6 

S-110.4 

DO  4  I*1,NL+1 
ZMSL ( I ) * ( 1-1 ) *DELZ 

HGEO < I ) -GAMMA* ( (RERTH*ZM8L (I ) ) / (RERTH+ZMBL (I) ) ) 

IF (HGEO ( I ) . LE. 1 1000. ) THEN 
XLMB— 6.5 
TB«288. 15 
HB*0. 

PB-101325. 

ELSEIF (HGEO ( I ) . GT. 1 1000. ) THEN 
XLMB-O. 

TB-216.65 
HB* 11000. 

PB-22632. 

ENDIF 

TZ ( I ) -TB+XLMB* ( HGEO ( I > -HB ) / 1 000 . 

IF (XLMB. NE.O.) THEN 

PR ( I ) -PB* (TB/TZ (I) ) ** (GOP*XMO* 1000. / (RSTAR*XLMB) ) 

ELSE 

PR  ( I ) -PB*EXP ( (-GOP*XMO* (HGEO ( I ) -HB) ) / (RSTAR*TB) ) 

ENDIF 

DENS ( I ) - (PR ( I ) *XMO) / (RSTAR*TZ  (I) ) 

VISK ( I ) -BETA* (TZ ( I ) ** 1 . 5) / <TZ  < I ) +S> /DENS ( I ) 

4  CONTINUE 


oono  n  nnnn  nnnn 


DO  9  M*1,NL+1 

DENS  <M> -DENS  (MX.  001 

VISK<M)«VISK(M) *10000. 

9  CONTINUE 

DAV I ES-MCD0NALD  METHOD  TO  COMPUTE  TERMINAL  VELOCITIES 
WT  »>  CDIR2  >»  RE  >»  VZ 


DO  13  1*1, NP 
SHRX-0. 

SHRY-0. 

COLRAD-CC 

RADLON-RR 

WT(I)*(4. /3. ) *3. 14159* ( (PDIA(I) /2. >  **3. ) *PDEN*980. 

DATA  TR/NL*0. / 

XM-0. 

YM-0. 

XR-RADLON 

YR-COLRAD 

IH  DESIGNATES  A  TRACE  PARTICLE’S  STARTING  LEVEL  NO.  IN  ATMOSPHERE 
LEVEL  1  *  SEA  LEVEL 
(IH-1)  X  DELZ  i  STARTING  HEIGHT 
IH* (NL+1) * (ZPDIA (I > /ZMAX) 

A*ZPDIA(I)-(IH-1)*DELZ 
IF (A.GE.DELZ/2. )  IH-IH+1 

TFALL-O. 

DO  12  J*IH, 1,-1 

COMPUTE  RESIDENCE  TIME  IN  EACH  LAYER 
(BETWEEN  TWO  LEVELS) 

0*8. *WT ( I) / (3. 1416*DENS ( J) * (VISK ( J) **2. ) > 

IF (Q. LT. 140. ) THEN 

R*Q/24. - (2. 3363E-4) *0**2. + (2. 0154E-6) *0**3. 

C-(6.9105E-9) *0**4. 

EL8EIF(Q.GE. 140.) THEN 

RL— 1 . 29536+ .986*(L0G10(Q))-.046677*((L0G10(Q)**2.))+ 

C. 0011235* ( (LOGIO(Q) **3. ) ) 

R*(10. )**RL 
END  IF 

V1*R»VISK ( J) /PDIA ( I ) 

IF (J.EQ. IH)THEN 
V2*V1 
60  TO  12 
END  IF 

VZA*(V2+Vl)/2. 

DENA* (DEN8(J) +DENS ( J+l ) ) /2. 


on  n  n  onon  non 


BLIPF-1 .  +  (2. 33E-7) / (PDIA  < I ) *. 01 IDENA* 1000. ) 

VZA-VZAI8LIPF 

V2-V1 

C  TR(J)  IS  RESIDENCE  TIME  IN  A  LAYER  (SECONDS) 
TR(J) -(DELZI100. ) /VZA 
TFALL«TFALL+TR(J) 


THERE  ARE  IH-1  LAYERS  OF  ATMOSPHERE  BELOM  A  PARTICLE 
12  CONTINUE 

DO  14  JJ-IH-1,1,-1 
US-0. 

VS-0. 

COMPUTE  DISTANCE  IN  BOTH  RADIANS  AND  METERS 
ZZ  IS  STARTIN6  HEIGHT  OF  PARTICLE 

ZZ-JJtDELZ 

USE  ZZ  TO  FIND  LEVELS  OF  WIND  COEFFICIENTS  ABOVE  AND  BELOW  ZZ 
CALL  LAYERS  (ZZ, LAYER, HL,HU) 

IF (LAYER. EQ. 1) THEN 
LVL-1 

CALL  UVCOMP (COLRAD, RADLON, LVL, US, VS, EPS, CU, CV) 

UU-US 

VU-VS 

UL-O. 

VL-O. 

GO  TO  19 

ELSEIF (LAYER. EQ. 13) THEN 
LVL-12 

CALL  UVCOMP (COLRAD, RADLON, LVL, US, VS, EPS, CU, CV) 

UU-US 
VU-VS 
UL-US 
VL-VS 
GO  TO  19 
ENDIF 


LVL-LAYER 

CALL  UVCOMP (COLRAD, RADLON, LVL, US, VS, EPS, CU, CV) 

UU-US 

VU-VS 

LVL-LVL-1 

CALL  UVCOMP (COLRAD, RADLON, LVL, US, VS, EPS, CU, CV) 

UL-US 

VL-VS 

INTERPOLATE  BETWEEN  THE  TWO  SPECTRAL  WIND  VECTORS. 
19  AA- (ZZ-HL) / (HU-HL) 

BB- (HU-ZZ) / (HU-HL) 

US-AAIUU+BBIUL 

VS-AAtVU+BBIVL 


IFtJJ.EQ.  IH-DTHEN 

USOLD-US 

VSOLD-VS 

60  TO  20 

ELSE 

SHRX»SHRX+<  ( (USOLD-US) /DELZ)tl2  )*TR(JJ+1) 

SHRY-SHRY+<  < (VSOLD-VS) /DELZ> *12  >«TR(JJ+1> 

USOLD-US 

VSOLD-VS 

ENDIF 


20  DYM-VSITR(JJ) 

DYR-DYh/RERTH 

YM-YM+DYM 

YR-YR+DYR 

VY  IS  POSITIVE  NORTHWARD 
COLRAD  IS  POSITIVE  SOUTHWARD 
COLRAD-COLRAD-DYR 

DXM-UStTR(JJ) 

DXR-DXM/ ( RERTHt  DS I N  <  COLRAD ) ) 

XM-XM+DXM 

XR-XR+DXR 

RADLON-RADLON+DXR 


IFtJJ.EQ. 1)THEN 

SHRX-SQRT (SHRX/TFALL) 
SHRY-SQRT <SHRY/TFALL ) 


PRAD-PDIA ( I ) /2. f . 01 

WRITE  <21 , 13) PRAD, TFALL, XM, SHRX, YM, SHRY 
FORMAT  (6E12.3) 

ENDIF 

CONTINUE 

CONTINUE 
GOTO  29 
CONTINUE 


SUBROUTINE  LAYERS  (ZZ, LAYER, HL, HU) 

THIS  ROUTINE  SETS  THE  WIND  LAYER  FOR  SELECTING 
SPECTRAL  COEFFICIENTS!  HEIGHTS  ARE  IN  METERS. 
IF  <ZZ . GE. 0. AND. ZZ. LT. 100) THEN 
LAYER- 1 
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HU- 100. 

ELSEIF(ZZ.GE. 100. AND. ZZ.LT. 1450)THEN 

LAYER-2 

HL-100. 

HU- 1450. 

ELSE I F ( Z  Z . GE . 1 450 . AND . Z  Z . LT . 3000 ) THEN 

LAYER-3  A 

HL-1450. 

HU-3000. 

ELSE IF ( Z  Z . GE . 3000 . AND . Z  Z . LT . 5600 ) THEN 

LAYER-4 

HL-3000. 

HU-5600. 

ELSE I F ( Z  Z . GE . 5600 . AND . Z  Z . LT . 7200 ) THEN 

LAYER-5 

HL-5600. 

HU-7200. 

ELSEIF (ZZ. GE. 7200. AND. ZZ.LT. 9200) THEN 

LAYER-6 

HL-7200. 

HU-9200. 

ELSEIF (ZZ. GE. 9200. AND. ZZ.LT. 10400) THEN 

LAYER-7 

HL-9200. 

HU-10400. 

ELSEIF (ZZ. GE. 10400.  AND. ZZ.LT. 1 1800) THEN 

LAYER-8 

HL-10400. 

HU-11800. 

ELSEIF (ZZ. GE. 1 1800. AND. ZZ.LT. 13600) THEN 

LAYER-9 

HL-11800. 

HU-13600. 

ELSEIF (ZZ.GE. 13600. AND. ZZ.LT. 16200) THEN 

LAYER-10 

HL- 13600. 

HU-16200. 

ELSEIF (ZZ.GE. 16200. AND. ZZ.LT. 18500) THEN 

LAYER- 11 

HL-16200. 

HU-18500. 

ELSEIF (ZZ.GE. 18500. AND. ZZ.LT. 20600) THEN 

LAYER- 12 

HL-18500. 

HU-20600. 

ELSEIF ( ZZ . GE. 20600. ) THEN 

LAYER- 12 

HL-20600. 

HU-100000. 

ENDIF 

RETURN 


on  nonooo  nnonnonnnnnnonon  no 


tmimmmmmmttmmmmutmmtmttimi 

SUBROUTINE  UVCOMP  ( COLRAD , RADLON, L VL , US, VS, EPS , CU , CV ) 

THIS  ROUTINE  COMPUTES  WIND  VECTOR  COMPONENTS 
FROM  NWS  SPECTRAL  COEFFICIENTS. 

THERE  ARE  JCAP+1  ZONAL  WAVE  NUMBERS  AND  JCAP+2  ORDINAL 
WAVE  NUMBERS  IN  THE  SPHERICAL  HARMONICS  SUMMATION. 

CU  AND  CV  ARE  COMPLEX  SPECTRAL  COEFFICIENTS 
FOR  W-E  AND  S-N  WIND  COMPONENTS. 

THERE  ARE  12  SETS  OF  COEFFICIENTS: 

ONE  SET  FOR  EACH  LEVEL  OF  ATMOSPHERE. 

EACH  SET  CONTAINS  992  COEFFICIENTS  (COMPLEX  PAIRS) 

FOR  JCAP-30. 

PARAMETER  (JCAP-30) 

COMPLEX  CU ( 1 i 12, 1 i JCAP+2, 1 i JCAP+1 ) , CV ( 1 l 12, 1 i JCAP+2, 1 i JCAP+1 ) 

COMPLEX  EIL,USUM,V8UM 

DOUBLE  PRECISION  EPS (JCAP+2, JCAP+1) 

DOUBLE  PRECISION  RADLON 
DOUBLE  PRECISION  COLRAD 
REAL  PLN (JCAP+2, JCAP+1) 

PI-3. 14159265 

COLRAD  AND  RADLON  ARE  COLATITUDE  AND  LONGITUDE  IN  RADIANS 


JCAP1-JCAP+1 

JCAP2-JCAP+2 


CALL  PLN3 (PLN, COLRAD, JCAP, EPS) 

C 

USUM-0. 

VSUM-0. 

DO  10  LL-1, JCAP1 
L-LL-1 

EIL-CMPLX (DCOS (LtRADLON) , DSIN(LtRADLON) ) 

DO  10  NN« 1 , JCAP2 

AA-2. 

IF(L.EQ.O)  AA-1. 

U8UM-U8UM+ AA IPLN (NN, LL) *CU (LVL, NN, LL) *EIL 
VSUM-VSUM+AAIPLN ( NN , LL ) *CV (LVL , NN , LL  >  *E I L 
10  CONTINUE 
C 

IF (COLRAD. LT. l.D- 10)  THEN 


onooo  no  on  n  o  oonnn  ono 


PRINT*,'  ***«  NORTH  POLE  «***' 

SO  TO  9999 
END  IF 
C 

C  WIND  VECTOR  COMPONENTS 

USsREAL (USUM) /DSIN (COLRAD) 

VS-REAL (VSUM) /DSIN (COLRAD) 

C 

9999  RETURN 
END 

***************************************************** 

SUBROUTINE  EPSLN(EPS, JCAP) 

*****  EPSILON  COEFFICIENTS  FOR  COMPUTING 
ASSOCIATED  LEGENDRE  POLYNOMIALS  WITH 
RECURSION  RELATION  IN  NWS30  P.24 

DOUBLE  PRECISION  EPS(l) 

DOUBLE  PRECISION  A 

JCAP-30 
JCAP 1 -JCAP* 1 
JCAP2-JCAP+2 

DO  100  LL-1.JCAP1 
L*LL-1 
JLE-L*JCAP2 

EPS ( JLE+INDE) -DSQRT (A) 

100  CONTINUE 

DO  200  LL-1.JCAP1 
JLE* (LL-1 ) *JCAP2 
EPS(JLE+l)-0.D+0 
200  CONTINUE 
RETURN 
END 

***************************************************** 

SUBROUTINE  PLN3 (PLN , COLRAD , JCAP , EPS ) 

*****  ASSOCIATED  LEGENDRE  POLYNOMIALS  COMPUTED 
WITH  RECURSION  RELATIONS  IN  BELOUSOV 
AND  NWS30  P.24 

DOUBLE  PRECISION  COLRAD 
DOUBLE  PRECISION  EPS(l) 

DOUBLE  PRECISION  8INLAT 
DOUBLE  PRECISION  C0S2 
DOUBLE  PRECISION  PROD 
DOUBLE  PRECISION  A 


DOUBLE  PRECISION  B 
DOUBLE  PRECISION  FL 
DOUBLE  PRECISION  PI 
DOUBLE  PRECISION  P2 
DOUBLE  PRECISION  P3 
DOUBLE  PRECISION  UNFLOW 

REAL  PLN(l) 

DATA  UNFLOW/. 75D-73/ 

SINLAT-DCOS(COLRAD) 

C0S2« l.D+O-SINLATISI NLAT 

PR0D*l.D+0 

A-l.D+O 

b-o.d+o 

JCAP1-JCAP+1 

JCAP2-JCAP+2 

DO  300  LL-1, JCAP1 

L-LL-1 

FL«L 

JLE-LIJCAP2 
IF(L.EQ.O)  60  TO  400 
A»A+2.D+0 
B-B+2.D+0 

FIX  UNDERFLOW 

IF (PROD. LE. UNFLOW)  PROD-O. D+0 

PR0D«PR0D*C0S2*A/B 
400  CONTINUE 

Pl-DSQRT (0. 5D+0tPR0D) 

PLN ( JLE+1 ) *SNGL (PI ) 

P2-DS0RT (2. D+0tFL+3. D+0) *SINLAT*P1 
PLN ( JLE+2) "SNGL (P2) 

DO  300  N-3.JCAP2 
LINDEX-JLE+N 

P3- (8INLATIP2-EPS (LINDEX-1 )*P1> /EPS (LINDEX) 
PLN (LINDEX ) »SNGL (P3) 

P1-P2 

P2«P3 

500  CONTINUE 

300  CONTINUE 
RETURN 


Appendix  G 


Dose  Program  Listing 


The  dose  program  listing  follows.  All  variables  and  routines  are 
defined  in  the.  program  header  blocks. 


PROGRAM  D0SE5 

********************************************************************** 
*  THIS  PROGRAM  CALCULATES  THE  DOSE  DUE  TO  THE  HOTLINE  INFORMATION 
♦PROVIDED  BY  HOLINE 


VARIABLES 

YO  *  YIELD  OLD(MT)  THE  YIELD  FROM  THE  PREVIOUS  HOTLINE 

YM  =■  YIELD  IN  MEGATONS  OF  THE  CURRENT  HOTLINE 

X  =  THE  TRACE  PARTICLE  X  COORDINATE 

Y  -  THE  TRACE  PARTICLE  Y  COORDINATE 

TFALL  -  THE  TRACE  PARTICLE  TIME  OF  FALL 

SHBK  -  THE  WIND  SHEAR  IN  THE  X  DIRECTION 

SHRY  »  THE  WIND  SHEAR  IN  THE  Y  DIRECTION 

R  -  TRACE  PARTICLE  RADIUS 

NP  =  THE  NUMBER  OF  TRACE  PARTICLES  FOR  EACH  HOTLINE 
BLAT  =  THE  BURST  POINT  LATITUDE 
BLON  =  THE  BURST  POINT  LONGTITUDE 

MINLAT  =  THE  MINIMUM  LATITUDE  FOR  EACH  HOTLINE  THAT  ACTIVITY  IS 
DEPOSITED  AT 

MAXLAT  =  THE  MAXIMUM  LATITUDE 
MINLON  =  THE  MINIMUM  LONGTITUDE 
MAXLON  =  THE  MAXIMUM  LONGTITUDE 

HC  =  THE  STABILIZED  CLOUD  HEIGHT  ACCORDING  TO  WSEG 
C  =  THE  LAURENT  SERIES  COEFFICIENTS  FOR  R  AS  A  FUNCTION  OF  TIME 
SZ  *  THE  STABILIZED  CLOUD  STANDARD  DEVIATION  IN  THE  VERTICAL 
ACCORDING  TO  WSEG 

SO  =  THE  STABILIZED  CLOUD  STANDARD  DEVIATION  IN  THE  HORIZONTAL 
ACCORDIit  TO  WSEG 

CELL  -  THE  DOSE  INFORMATION  FOR  EACH  1  DEG  LAT  BY  1  DEG  LONG 
CELL  IN  THE  US 
DIMENSIONED  (I,J,K) 

I  -  LATITUDE 
J  -  LONGTITUDE 

X-l  -  DOSE  TO  INFINITY  (RADS-TISSUE) 

K-2  -  UNIT  TIME  REFERENCE  DOSE  RATE 
K-3  -  DOSE  RATE  AT  1  DAY 
£-4  -  DOSE  RATE  AT  7  DAYS 


K-5  -  DOSE  RATE  AT  14  DAYS 
K-6  -  DOSE  RATE  AT  30  DAYS 
CITY  *  THE  DOSE  TO  INFINITY  FOR  EACH  CITY 
NKILL  *  THE  NUMBER  OF  US  POPULATION  DEATHS  DUE  TO  FALLOUT 
A  *  THE  SLOPE  OF  EACH  HOTLINE  SEGMENT 
B  -  THE  INTERCEPT  OF  EACH  HOTLINE  SEGMENT 
NBURS  -  THE  NUMBER  OF  BURSTS 
HLINE, WLINE  >  THE  HEIGHT (NS)  AND  WIDTH (EW)  OF  AN  AREA  OF  BURSTS 
FV  -  VOLUME  FRACTION  OF  ACTIVITY 

IF  FV>1  THEN  THE  2.5  MOMENT  IS  USED  FOR  ACTIVITY  SIZE  DIST. 

FF  -  THE  FISSION  FRACTION  OF  THE  YIELD 

ALPHA  «  THE  LOG  MEAN  OF  THE  ACTIVITY  SIZE  DISTRIBUTION 

BETA  -  THE  LOG  SLOPE  OF  THE  ACTIVITY  SIZE  DISTRIBUTION 


SUBROUTINES  CALLED 

UNITS  -  CONVERTS  TFALL  FROM  SECONDS  TO  HOURS,  SHEAR  VALUES  FROM 
PER  SECOND  TO  PER  HOUR,  AND  <X,Y>  FROM  METERS  TO 
KILOMETERS 

VALUES  -  FINDS  THE  STABILIZED  CLOUD  YIELD  DEPENDENT  VARIABLES 
LIMITS  -  FINDS  THE  BOX  CONTAINING  THE  HOTLINE  AND  THE  CELLS 
AFFECTED  BY  THE  HOTLINE 

FALOUT  -  FINDS  THE  DOSE  INFORMATION  FOR  EACH  CELL  AND  CITY 

DEATH  -  FINDS  THE  POPULATION  DEATHS  DUE  TO  THE  FALLOUT 

RITE  -  OUTPUT  THE  NUMBER  OF  DEATHS  AND  THE  CELL  DOSE  INFORMATION 


FILES  USED 
INPUT 

TAPE21  -  THIS  FILE  IS  CREATED  BY  THE  HOTLINE  LOCATOR,  AND  FOR  EACH 
HOTLINE  HASi  THE  BURST  LOCATION,  YIELD,  NUMBER  OF  BURSTS 
SIZE  OF  THE  BURST  AREA,  VOLUME  AND  FISSION  FRACTIONS,  AND 
ACTIVITY  SIZE  DISTRIBUTION  ALPHA  AND  BETA.  AND  FOR  EACH 
TRACE  PARTICLE  HAS  THE  PARTICLE  TIME  OF  FALL,  SIZE,  THE 
WIND  SHEAR  AND  THE  <X,Y>  LOCATION 

OUTPUT 

TAPE22  -  THIS  IS  THE  PROGRAMS  OUTPUT  FILE  CONTAINING  THE  NUMBER  OF 
POPULATION  DEATHS  AND  THE  DOSE  INFORMATION  FOR  EACH  CELL 

PARAMETER (NP- 10) 

REAL  YO, YM, X(OiNP) ,Y (Oi NP) , TFALL (Oi NP) , SHRX (Oi NP) , 

C  8HRY (Oi NP) , R (Oi NP) , BLAT, BLON, M INLAT, MAXLAT , MINLON, MAXLON, 

C  HC,C(7) ,SZ,S0,CELL(25i49,67l 124,6) , CITY (316) ,NKILL,A(OiNP-l) , 

C  B (0 i NP- 1 ) , NBURS , HL INE , WL I NE, FV, FF , ALPHA , BETA 

INTEGER  I 

DATA  CELL/870040./ 

DATA  CITY/31640./ 

DATA  YO/O. / 

REWIND  21 
5  CONTINUE 


READ (21 , 1 1 , END-999) BLAT, BLON, YM, NBURS, HLINE, WLINE, FF, FV, 


C  ALPHA, BETA 
DO  10  I*NP, 1,-1 

READ (21 , 12) R ( I ) , TFALL ( I ) , X ( I ) , SHRX ( I ) , Y ( I ) , SHRY ( I > 

10  CONTINUE 

1 1  FORMAT  (F5. 2, F8. 2, El  1 . 3, 3F5. 0,  F5. 2,  F6. 2, 2F10. 5,  F5. 1 ) 

12  FORMAT (6E12. 3) 

CALL  UNITS (TFALL, X , SHRX , Y, SHRY , NP) 

CALL  LIMITS(X,Y,BLAT, BLON , M I NL AT , MA X LAT , M I NLON , 

C  MAXLON, NP, HLINE, WLINE) 

IF (YM.NE.YO)THEN 
CALL  VALUES (YM,HC,C,SZ, SO) 

YQ“YM 

ENDIF 

CALL  FALOUT (CELL, CITY, X, Y, R,  TFALL, YM, BLAT, BLON, MAXLAT, 
C  MINLAT, MAXLON, MINLON, SO, SZ , C, NP, SHRX , SHRY, A, B, 

C  NBURS, HLINE, WLINE, FF, FV, ALPHA, BETA) 

GOTO  3 

999  CONTINUE 

CALL  DEATH (CELL, CITY, NKILL) 

CALL  RITE (CELL, NKILL) 

END 


SUBROUTINE  UNITS (TFALL, X, SHRX, Y, SHRY, NP) 

mmmmmmmmiummmmmmtmmmmmmmu 

THIS  ROUTINE  CHANGES  THE  UNITS  OF  THE  VARIABLES  IN  THE  ARGUMENT 
LIST 


VARIABLES 

TFALL  -  TIME  OF  FALL,  SECONDS  TO  HOURS 
X  -  TRACE  PARTICLE  X  COORDINATE,  METERS  TO  KILOMETERS 
SHRX  -  WIND  SHEAR  IN  X  DIRECTION,  /SEC  TO  /HOUR 
Y  «  TRACE  PARTICLE  Y  COORDINATE,  METERS  TO  KILOMETERS 
SHRY  -  WIND  SHEAR  IN  Y  DIRECTION,  /SEC  TO  /HOUR 
NP  -  NUMBER  OF  TRACE  PARTICLES 

...NOTE... THE  IF  STATEMENT  INSURES  A  FINITE  SLOPE  FOR  EACH 
HOTLINE  SEGMENT 


NO  SUBROUTINES  CALLED 


DO  10  1*1, NP 

TFALL < I ) -TFALL < I ) /3600. 

X(I)-X(I)/1000. 

Y(I)*Y(I)/1000. 

SHRX(I)-SHRX(I)*3600. 

SHRY ( I ) *SHRY ( I >  *3600. 

10  CONTINUE 

X(0)»0. 

Y (0)*0. 

TFALL (0)*0. 

SHRX (0)*0. 

SHRY(0)=0. 

DO  20  I*0,NP-1 

IF (ABS(X(I+1)-X  (I) )  .LE.  l.E-03)  X(I+1)*X  (I+D+.OOl 
20  CONTINUE 

END 

SUBROUTINE  VALUES <YM,HC,C,SZ, SO) 

tmummmmtmmmummmmmmmmttmmitmm 

t  THIS  ROUTINE  COMPUTES  YIELD  DEPENDENT  VALUES  ACCORDINS  TO  WSEG 
*  FORMULAS 


VARIABLES 

YM  -  YIELD  IN  MEGATONS 

HC  -  STABILIZED  CLOUD  HEIGHT  IN  KILOMETERS 
C  -  LAURENT  COEFFICIENTS  FOR  R(T) 

CC  *  LAURENT  COEFFICIENTS  FOR  NEXT  HIGHER  ALTITUDE 
SZ  ■  STABILIZED  CLOUD  STANDARD  DEVIATION  IN  THE  VERTICAL  IN 
KILOMETERS 

SO  *  STABILIZED  CLOUD  STANDARD  DEVIATION  IN  THE  HORIZONTAL  IN 
KILOMETERS 

Z1,Z2  *  HEIGHTS  FROM  THE  LAURENT  COEFFICIENT  TABLE  FOR 
INTERPOLATION 

YL  *  LOG  OF  YIELD  IN  MEGATONS 

FRAC  *  LINEAR  INTERPOLATION  FACTOR  BETWEEN  ALTITUDES 


NO  SUBROUTINES  ARE  CALLED 


FILES 
INPUT 

TAPE30  -  A  TABLE  OF  COLARCO'S  LAURENT  COEFFICIENTS  AS  A  FUNCTION 
OF  STABILIZED  CLOUD  HEIGHT,  TABULATED  EVERY  200  METERS 

********************************************************************* 

REAL  YM, HC, C (7) , CC  <7> , SZ, SO,  Z 1 ,  Z2,  YL 
INTEGER  I 


YL*LOG(YM) 

HC* (44. +6. 1*YL-. 205* (YL+2. 42) tABS (YL+2. 42) >  *1 . 609/5.280 
8Z*. 18*HC 


S0*1 . 609IEXP < . 70+1 . /3. *YL-3. 25/ (4.0+ (YL+5. 41*2) ) ) 


REWIND  30 

READ (30,100)21, (C< I), 1*1,4) 

READ (30, 100)21, <C(I),I=5,7) 

10  CONTINUE 

READ (30, 100) 22, (CC(I) , 1*1 , 4) 

READ (30, 100) 22, (CC(I) , 1*5,7) 

100  FORMAT (F5.1,4E1 1.5) 

IF < 22. LT. HC) THEN 
21*22 

DO  20  1*1,7 
C(I)*CC(I) 

20  CONTINUE 
GOTO  10 
ENDIF 

FRAC*(HC-21)/ (22-21) 

DO  30  1-1,7 

C ( I ) *FRACt  <CC ( I ) — C ( I ) ) +C ( I ) 

30  CONTINUE 

END 

SUBROUTINE  LIMITS ( X , Y, BLAT, BLQN, MINLAT, MAXLAT, 

C  MINLON,MAXLON,NP,HLINE, WLINE) 

mmmmmtmmmmmmmmmmmtmmmmtttttm 

*  THIS  ROUTINE  DEFINES  THE  LIMITS  OF  THE  LATITUDE  LONGTITUDE  BOX 

*  WHICH  CONTAINS  THE  HOTLINE  AND  ALL  OF  THE  CELLS  AFFECTED  BY  THE 

*  HOTLINE 


VARIABLES 

X  *  TRACE  PARTICLE  X  COORDINATE  IN  KILOMETERS 

Y  -  TRACE  PARTICLE  Y  COORDINATE  IN  KILOMETERS 

BLAT  *  BURST  LATITUDE 

BLON  *  BURST  LONGTITUDE 

MAXLAT  *  MAXIMUM  LATITUDE  OF  THE  BOX 

MINLAT  -  MINIMUM  LATITUDE  OF  THE  BOX 

MINLON  «  MINIMUM  LONGTITUDE  OF  THE  BOX 

MAXLON  *  MAXIMUM  LONGTITUDE  OF  THE  BOX 

NP  -  THE  NUMBER  OF  TRACE  PARTICLES 

MAXX  -  MAXIMUM  X  OF  BOX 

MINX  -  MINIMUM  X  OF  BOX 

MAXY  -  MAXIMUM  Y  OF  BOX 

MINY  -  MINIMUM  Y  OF  BOX 

HLINE  -  HEIGHT  OF  AREA  OF  BURSTS 

WLINE  -  WIDTH  OF  AREA  OF  BURSTS 


NO  SUBROUTINES  ARE  CALLED 


NO  FILES  ARE  USED 

mitiitmmmtitimmmmmmttmmimmtmmtmim 
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.  V  •  *  . 


PARAMETER  <  P I =3 . 1 4 1 59265 ) 

REAL  X (0: NP) , Y (Os  NP) , BLAT, BLON, MINLAT, MAXLAT, 

C  MINLON, MAXLON, MAX  X , MINX , MAX Y, MINY, HL INE, WL INE 

INTEGER  I 

MAXX*0. 

MINX=0. 

MAXY=0. 

MINY=0. 

DO  10  1=0, NP 

MAXX=MAX (X (I) ,MAXX) 

MINX=MIN (X (I ) , MINX) 

MAXY*MAX (Y ( I > , MAXY> 

MINY=MIN(Y (I) ,MINY) 

10  CONTINUE 

MAXY-MAXY+HL1NE/2. 

MINY-MINY-HLINE/2. 

MAXX=MAXX+WLINE/2. 

MINX-MINX-WLINE/2. 


MAXLAT=MAXY/ 1 10. 94  +  BLAT  +6. 

MINLAT=MINY/1 10. 94  +  BLAT  -6. 

MAXLON=-MINX/ (COS (BLATIPI/1B0. ) tl 10. 94) +  BLON  +B. 

MINLON=-MAXX/ (COS (BLAT*PI/180.) *110.94)+  BLON  -B. 

MAXLAT=MIN (MAXLAT, 50. ) 

MINLAT =MAX (MINLAT, 25. ) 

MAXLON*MIN (MAXLON, 125. ) 

MINLON=MAX (MINLON, 67. ) 

END 

SUBROUT  I NE  FALOUT  ( CELL ,  C I TY ,  X ,  Y ,  R ,  TF  ALL ,  YM ,  BLAT ,  BLON ,  MAXLAT , 

C  MINLAT, MAXLON, MINLON, SO,  SZ , C, NP , SHRX , SHRY, A, B, 

C  NBURS , HL I NE ,  WL I NE ,  FF , F V , ALPHA , BET A ) 

ll*tll**t***ll*l***l****l***t*t ***********************  *************** 
THIS  ROUTINE  CALCULATES  THE  DOSE  FOR  EACH  CELL  AND  CITY  AFFECTED 
BY  THE  HOTLINE. 


VARIABLES 

CELL  ■  THE  DOSE  INFORMATION  FOR  EACH  CELL 
DIMENSIONED  <I,J,K) 

I  -  CELL  LATITUDE 
J  -  CELL  LONGTITUDE 
K«1  DOSE  TO  INFINITY 
K»2  UNIT  TIME  REFERENCE  DOSE  RATE 
K«3  DOSE  RATE  AT  1  DAY 
K-4  DOSE  RATE  AT  7  DAYS 
K=5  DOSE  RATE  AT  14  DAYS 
K-6  DOSE  RATE  AT  30  DAYS 
CITY  -  DOSE  TO  INFINITY  FOR  EACH  CITY 
X  =  TRACE  PARTICLE  X  COORDINATE  IN  KILOMETERS 


Y  -  TRACE  PARTICLE  Y  COORDINATE  IN  KILOMETERS 

R  «  TRACE  PARTICLE  RADIUS  IN  METERS 

TFALL  -  TRACE  PARTICLE  TIME  OF  FALL  IN  HOURS 

YM  ■  YIELD  IN  MEGATONS 

BLAT  «  BURST  LATITUDE 

BLON  -  BURST  LONGTITUDE 

MAXLAT  -  BOX  MAXIMUM  LATITUDE 

MINLAT  -  BOX  MINIMUM  LATITUDE 

MAXLON  -  BOX  MAXIMUM  LONGTITDE 

MINLON  «  BOX  MINIMUM  LONGTITDE 

SO  >  STABILIZED  CLOUD  HORIZONTAL  STANDARD  DEVIATION 
SZ  -  STABILIZED  CLOUD  VERTICAL  STANDARD  DEVIATION 
C  «  LAURENT  COEFFICIENTS  FOR  R  AS  A  FUNCTION  OF  TIME 
SHRX  -  WIND  SHEAR  IN  X  DIRECTION  FOR  EACH  TRACE  PARTICLE 
SHRY  -  WIND  SHEAR  IN  Y  DIRECTION  FOR  EACH  TRACE  PARTICLE 
IMP  «  NUMBER  OF  TRACE  PARTICLES 

HOTX,HOTY  «  POSITION  ON  THE  HOTLINE  DEPOSITING  ACTIVITY  AT 
(XX, YY) 

HOTTA  «  THE  TIME  OF  ARRIVAL  AT  (HOTX.HOTY) 

HOTSX.HOTSY  *  THE  WIND  SHEAR  AT  (HOTX, HOTY) 

DINF  *  TIME  INTEGRATED  DOSE  TO  INFINITY 

UTRDR  *  UNIT  TIME  REFERENCE  DOSE  RATE 

DR1  *  DOSE  RATE  AT  1  DAY 

DR7  *  DOSE  RATE  AT  7  DAYS 

DR14  -  DOSE  RATE  AT  14  DAYS 

DR30  «  DOSE  RATE  AT  30  DAYS 

CL AT  =  CITY  LATITUDE 

CLON  -  CITY  LONGTITUDE 

LON  -  CELL  CENTER  LONGTITUDE 

LAT  -  CELL  CENTER  LATITUDE 

CLOMIN  =  MINIMUM  CELL  CENTER  LONGTITUDE  OF  FALLOUT  BOX 
CLOMAX  =■  MAXIMUM  CELL  CENTER  LONGTITUDE 
CLAM IN  «  THE  MINIMUM  CELL  CENTER  LATITUDE 
CLAMAX  *  MAXIMUM  CELL  CENTER  LATITUDE 

A,B  *  SLOPE  AND  INTERCEPT  OF  THE  STRAIGHT  LINE  EQUATION  FOR  EACH 
HOTLINE  SEGMENT 

SWI  -  LOGICAL  SWITCH;  IF  SWI*.TRUE.  THEN  THERE  IS  NO  POINT  ON 
THE  HOTLINE  WITH  A  CROSSWIND  VECTOR  TO  THE  CELL  CENTER  OR 
CITY  LOCATION.  NO  ACTIVITY  DEPOSITED  AT  THAT  POINT. 

CSWI  -  LOGICAL  SWITCH;  IF  CSWI-.TRUE.  THEN  FINDING  ACTIVITY  AT 
A  CITY 

XX, YY  -  POINT  AT  WHICH  FINDING  ACTIVITY  DEPOSITED 
NBURS  -  NUMBER  OF  BURSTS 

HLINE, WLINE  -  HEIGHT  AND  WIDTH  OF  AN  AREA  OF  BURSTS 
FV  «  VOLUME  FRACTION  OF  ACTIVITY 
FF  ■  FISSION  FRACTION  OF  YIELD 

ALPHA, BETA  ■  PARAMETERS  OF  ACTIVITY  SIZE  DISTRIBUTION 


SUBROUTINES  CALLED 

8L0INT  -  FINDS  THE  SLOPE  AND  INTERCEPT  FOR  THE  LINE  EQUATION  OF 
EACH  HOTLINE  SEGMENT 

POSIT  -  FINDS  THE  POSITION  (XX, YY)  FROM  THE  BURST  POINT  FOR  THE 
LAT  LON  COORDINATES 
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HLLOC  -  FINDS  THE  HOTLINE  COORDINATE  (HOTX, HOTY)  WHICH  DEFINES  THE 
ACTIVITY  DEPOSITED  AT  (XX, YY) 

SMEAR  -  FINDS  THE  DOSE  INFORMATION  AT  (XX, YY)  DUE  TO  <HOTX,HOTY> 


FILES  USED 
INPUT 

TAPE53  -  THIS  FILE  IS  A  LIST  OF  EACH  CITY'S  LATITUDE,  LONGTITUDE 
AND  POPULATION 

mmmsttmmmmttmmmmmmmmmmmmmmt 

REAL  CELL (25: 49, 67: 124,6) , CITY <31 6) , X (0:NP) , Y(0:NP) ,R(0:NP) , 

C  TFALL (0: NP) , YM,  BLAT, BLON, MAXLAT, MINLAT, MAXLON, MINLON, SO, SZ, 

C  XX, YY, LAT, LON, SHRX (0: NP) , SHRY (0: NP) , HOTX, HOTY, HOTTA, DINF, 

C  DR1,DR7,DR14,UTRDR,CLAT,CLQN,C(7),DR30,A(0:NP-1),B(0:NP-1>, 

C  NBURS,HLINE,WLINE, ALPHA, BETA 
INTEGER  I 
LOGICAL  SWI.CSWI 

SW Is. FALSE. 

CSWI-. FALSE. 

CALL  SLOINT (X,Y, A,B,NP) 

CLOMAX-INT (MAXLON) +. 5 
IF (CLOMAX . GT . 1 25 . ) CLOMAX- 124.5 
CLOMIN-INT (MINLON) +.5 
CLAMAX- I NT (MAXLAT) +.  5 
I F (CLAMAX . GT . 50 . ) CLAMAX  *49 . 5 
CLAMIN-INT (MINLAT) +.5 


DO  10  LON=CLOMIN, CLOMAX 
DO  10  LAT-CLAMIN, CLAMAX 
CALL  POS I T ( X  X ,  Y Y , BLAT ,  BLON , LAT , LON ) 

CALL  HLLOC  <  X  X ,  YY ,  X ,  Y , TFALL , HOTX , HOTY , HOTTA , SWI , SHRX , SHRY , 

C  HOTSX , HOTSY , NP, A, B) 

IF (SWI) THEN 
SWI*. FALSE. 

GOTO  10 
ENDIF 

CALL  SMEAR  (HOTX,  HOTY,  XX,  YY,  HOTTA,  DINF,  UTRDR,  DR1 ,  DR7,  DR14,  YM, 
C  SO , SZ , C , HOTSX , HOTSY , DR30 , NBURS , HL I NE , WL I NE , ALPHA , BETA , 

C  FV,FF,CSWI,BLAT) 

CELL ( INT (LAT) , INT (LON) , 1 ) -CELL ( INT (LAT) , INT (LON) , 1 ) +DINF 
CELL ( INT (LAT) , INT (LON) , 2) -CELL ( INT (LAT) , INT (LON) , 2) +UTRDR 
CELL ( INT (LAT) , INT (LON) , 3) -CELL  < INT (LAT) , INT (LON) , 3) +DR1 
CELL (INT (LAT) , INT (LON) , 4) -CELL ( INT (LAT) , INT (LON) , 4) +DR7 
CELL ( INT (LAT) ,  INT (LON) , 5) -CELL ( INT (LAT) , INT (LON) , 5) +DR14 
CELL  < INT (LAT) , INT (LON) , 6) -CELL (INT (LAT) , INT (LON) ,6) +DR30 
10  CONTINUE 


REWIND  53 
SWI-. FALSE 
CBWI-.TRUE 


DO  20  1-1,316 

READ (53, 1 1 )CLAT , CLON, POP 

IF (CLAT . 6T . MAXLAT+1  .OR.  CLAT . LT. MINLAT-1 .  .OR.  CLON.GT.MAXLON+1 
C  .OR.  CLON . LT . M I NLON- 1 . ) GOTO  20 

CALL  POSIT ( XX , YY , BLAT , BLON, CLAT , CLON) 

CALL  HLLOC ( X  X , YY , X , Y , TFALL , HOTX , HOTY , HOTTA , SW I , SHRX , SHRY , 

C  HOTSX , HOTSY ,  NP ,  A , B ) 

IF(SWI)THEN 
SWI-. FALSE. 

GOTO  20 
ENDIF 

CALL  SMEAR <HOTX, HOTY, XX, YY, HOTTA, DINF, UTRDR, DR1 , DR7, DR14, YM, 

C  SO , SZ , C , HOTSX , HOTS Y , DR30 , NBURS , HL I NE , ML I NE , ALPHA , 

C  BETA, FV,FF,CSWI, BLAT) 

CITY  < I ) -CITY ( I ) +DINF 
20  CONTINUE 
11  FORMAT (2F9. 2, FI 0.0) 

END 


SUBROUTINE  DEATH (CELL, CITY, NKILL) 

umimmtmmmmmttmmtmmmmmtmmummu 

i  THIS  ROUTINE  CALCULATES  THE  NUMBER  OF  DEATHS  DUE  TO  FALLOUT 


VARIABLES 

CELL  -  THE  DOSE  INFORMATION  FOR  EACH  1  DEG  LAT  BY  1  DEG  LONG 
CELL  IN  THE  US 
DIMENSIONED ( I, J,K) 

I  -  LATITUDE 
J  -  LONGTITUDE 
K-l  -  DOSE  TO  INFINITY 
CITY  -  THE  DOSE  TO  INFINITY  FOR  EACH  CITY 
NKILL  -  THE  TOTAL  NUMBER  OF  DEATHS 

DOSS  «  DOSE  SURE  SAFE,  NO  DEATHS  FOR  DOSES  LESS  THAN  THIS  VALUE 
DOSK  «  DOSE  SURE  KILL,  100'/.  DEATHS  FOR  DOSE  GREATER  THAN  THIS 
A.B  *  DUMMY  VARIABLES  FOR  THE  CITY  LAT  LON  WHICH  IS  NOT  USED 
POP  -  POPULATION  FOR  EACH  CELL  AND  THOUSANDS  FOR  EACH  CITY 
CNF  «  THE  PROBABLITY  OF  DEATH 
PF  «  PROTECTION  FACTOR  FROM  DOSE  GIVEN  BY  FALLOUT  SHELTER 
ALPHA  ■  LOG (MEAN)  OF  PROBABILITY  OF  DEATH  DISTRIBUTION 
BETA  -  LOG (SLOPE)  OF  PROBABILITY  OF  DEATH  DISTRIBUTION 
Z  -  ARGUMENT  FOR  A  CUMULATIVE  LOG-NORMAL  DISTRIBUTION 


SUBROUTINES  CALLED 

FUNCTION  CNF  -  CALCULATES  THE  PROBABILITY  OF  DEATH,  ASSUMING 
THE  PROBABILITY  IS  A  CUMULATIVE  LOG  NORMAL  FUNCTION,  OF 
DOSE,  AND  THAT  DOSS  HAS  A  PROBABLITY  OF  .01  AND  DOSK  HAS  A 
PROBABILTY  OF  .99 


FILE8  USED 

TAPE55  -  THE  RURAL  POPULATION  DATA  FOR  EACH  CELL 
TAPE53  -  THE  URBAN  POPULATION  FOR  EACH  CITY 
mmmmmmmtmmmMmmmtmmtmmmtmmm 


REAL  CELL (25s 49, 67i  124,6)  ,CITY (316) ,NKILL,DOSS,DOSK,PF,A,B, 
C  POP, ALPHA, BETA, Z 
INTEGER  I,J 

DATA  DOSS/ 169./ 

DATA  DOSK/1198.  / 

DATA  PF/. 333333/ 

ALPHA-. 5«L06 (DOSK* DOSS) 

BETA-LOG (DOSK/DOSS) /4. 656 

REWIND  55 
NK ILL-O. 

DO  20  1=25,49 
DO  20  J=67, 124 

READ  (55, 21 )  A,  B,  POP 
IF (CELL ( I , J , 1 ) *PF . LT . DOSS ) GOTO  20 
IF (CELL ( I , J,  1 ) *PF.  GT.  DOSK) THEN 
NK I LL=NK I LL+POP 
ELSE 

Z=  (LOG (CELL (I, J, 1) *PF) -ALPHA) /BETA 
NK I LL=NK I LL+PGP*CNF ( Z ) 

END  IF 
CONTINUE 

FORMAT (2F8.1.E10. 3) 

REWIND  53 
DO  30  1=1,316 
READ (53, 40) A, B, POP 
IF (CITY ( I ) *PF.  LT. DOSS) GOTO  30 
IF  (CITY ( I ) *PF. GT. DOSK) THEN 
NK I LL-NK I LL+POP* 1 000 . 

ELSE 

Z= (LOG (CITY ( I ) *PF) -ALPHA) /BETA 
NK I LL-NKI LL+POP* 1 000 . *CNF ( Z ) 

ENDIF 

CONTINUE 

FORMAT (2F9. 2, F10.0) 


SUBROUTINE  RITE (CELL, NKILL) 

********************************************************************** 

*  THIS  ROUTINE  WRITES  THE  NUMBER  OF  DEATHS  AND  THE  DOSE  INFORMATION 

*  FOR  EACH  CELL  TO  TAPE22 


VARIABLES 

CELL  =  THE  DOSE  INFORMATION  FOR  EACH  CELL 
DIMENSIONED  (I,J,K> 

1=1  -  CELL  LATITUDE 
J  -  CELL  LONGTITUDE 
K=1  DOSE  TO  INFINITY 


K-2  UNIT  TIME  REFERENCE  DOSE  RATE 
K-3  DOSE  RATE  AT  1  DA 
K-4  DOSE  RATE  AT  7  DAYS 
K=5  DOSE  RATE  AT  14  DAYS 
K*6  DOSE  RATE  AT  30  DAYS 
NKILL  *  THE  TOTAL  NUMBER  OF  DEATHS 
ND  *  CHARACTER  STRING  ’THE  NUMBER  OF  DEATHS  *  ’ 

HEADER  =  CHARACTER  STRING  HEADER  FOR  CELL  DOSE  INFORMATION 


NO  SUBROUTINES  ARE  CALLED 


*  FILES 
t  OUTPUT 

*  TAPE22  -  CONTAINS  THE  NUMBER  OF  DEATHS  DUE  TO  FALLOUT,  AND  THE 

«  DOSE  INFORMATION  FOR  EACH  CELL  CENTER 

mmmmmmmmmmtmmmmmtmmmmtmmmt 

REAL  CELL <25i 49,67i 124,6) , NKILL 
INTEGER  I,J 

CHARACTER  NDt22, HEADER475 
REWIND  22 

ND=’THE  NUMBER  OF  DEATHS  *’ 

HEADER*’  LAT  LONG  DINF  UTRDR  DR1  DR7' 

C  //  ’  DR14  DR30’ 

WRITE <22, 10) ND, NKILL 
WRITE  <22, 1DHEADER 
DO  20  1*25,49 
DO  20  J*67, 124 

NRITE<22,12)REAL(I)+.5,REAL<J)+.5, <CELL<I, J,K) ,K-1,6) 

20  CONTINUE 

10  FORMAT <A22,EB. 3) 

11  FORMAT <A75) 

12  FORMAT  <2(F5. 1,2X) , 6(E8.3,3X) ) 

END 

SUBROUT I NE  POS I T  <  X  X , YY , BL AT , BLON , L AT , LON ) 
iimmmmmtiimimmimiiimmmmmmmmmmm 
t  GIVEN  THE  BURST  LATITUDE  AND  THE  BURST  LONGTITUDE,  THIS  ROUTINE 
4  FINDS  THE  POSTION  <XX,YY)  CORRESPONDING  TO  <LAT,LON) 


VARIABLES 

XX  -  THE  DISTANCE  ALONG  THE  X  AXIS  FOR  POINT <LAT, LON) 
YY  -  THE  DISTANCE  ALONG  THE  Y  AXIS  FOR  POINT (LAT, LON) 
BLAT  -  BURST  LATITUDE  IN  DEGREES,  0  ON  THE  X  AXIS 
BLON  *  BURST  LONGTITUDE  IN  DEGREES,  0  ON  THE  Y  AXIS 
LAT  -  LATITUDE  OF  A  POINT 
LON  -  LONGTITUDE  OF  A  POINT 
NOTE i  THERE  ARE  ABOUT  110.94  KM/DEGREE  OF  LATITUDE 


*  NO  FILES  ARE  USED 


PARAMETER (PI-3. 14159265) 

REAL  XX , YY, BLAT, BLON , LAT , LON 

XX- (BLON-LON ) tCOS ( BL AT»P I / 1 BO . ) *  1 1 0 . 94 
YY- ( LAT-BLAT ) * 1 1 0 . 94 

END 

SUBROUT I NE  HLLOC ( X  X , Y Y ,  X ,  Y , TF ALL , HOT  X , HOT Y , HOTTA ,  SW I , SHRX , SHR Y , 
C  HOTSX,HOTSY,NP,A,B> 

THIS  ROUTINE  LOCATES  THE  POSITION  ON  THE  HOTLINE  WHICH  DETERMINES 
THE  UNIT  TIME  REFERENCE  DOSE  RATE  AT  THE  POSITION (XX, YY)  OFF  THE 
HOTLINE,  ACCORDING  TO  HOPKINS'  SMEAR  EQUATION.  THE  TIME  OF  ARRIVAL 
AND  WIND  SHEAR  AT  THE  HOTLINE  POINT  ARE  ALSO  FOUND. 


VARIABLES 

XX  -  X  COORDINATE  OF  A  POINT  OFF  THE  HOTLINE 
YY  -  Y  COORDINATE  OF  A  POINT  OFF  THE  HOTLINE 
X,Y  -  TRACE  PARTICLE  COORDINATES 
TFALL  -  TRACE  PARTICLE  TIMES  OF  FALL 

HOTX,HOTY  «  THE  HOTLINE  COORDINATE  TO  CALCULATE  ACTIVITY  AT 
(XX  YY) 

HOTTA  -  TIME  OF*  ARRIVAL  AT  (HOTX.HOTY) 

SWI  -  LOGICAL  VARIABLE,  WHEN  TRUE,  THERE  IS  NO  POINT  ON  THE 
HOTLINE  CONTRIBUTING  ACTIVITY  AT  (XX, YY) 

SHRX , SHRY  «  WIND  SHEAR  FOR  TRACE  PARTICLES  IN  X  AND  Y  DIRECTION 

HOTSX.HOTSY  -  WIND  SHEAR  AT  (HOTX,HOTY)  IN  X  AND  Y  DIRECTION 

NP  «  NUMBER  OF  TRACE  PARTICLES 

A  -  THE  SLOPE  OF  A  HOTLINE  SEGMENT 

B  -  THE  INTERCEPT  OF  A  HOTLINE  SEGMENT 

AA,BB,CC  -  QUADRATIC  COEFFICIENTS  TO  FIND  (HOTX.HOTY) 

XXI, XX2  -  TRIAL  HOTLINE  X  COORDINATES 
T  -  THE  PARAMETER  IN  THE  PARAMETRIC  REPRESENTATION  OF  THE 
HOTLINE  SEGMENT 


SUBROUTINES  CALLED 

LINSEG  -  DETERMINES  IF  A  TRIAL  X  IS  ON  A  HOTLINE  SEGMENT,  AND  IF 
BOTH  TRIAL  X  ARE  ON  A  HOTLINE  SEGMENT,  PICKS  THE 
EARLIEST  ONE 


NO  FILES  ARE  USED 

REAL  XX, YY, X <0« NP) , Y (0» NP) , TFALL <0» NP) , HOT X, HOT Y, HOTTA, 

C  SHRX (Oi NP) , SHRY (Oi  NP) , HOTSX.HOTSY, AA, BB, CC, XXI , XX2, 
C  T,A(OiNP-l),B(OiNP-l) 

LOGICAL  SWI 


INTEGER  I 


DO  10  I-O.NP-l 
AA— A(I)*«2-1. 

BB-XX+A ( I ) IYY-2. «A ( I ) *B ( I ) 

CC«B  < I >  tYY-B  < I ) $$2 
IF  (BBM2.LT.  4.  *AA*CC>  GOTO  10 
XX1*(-BB+SQRT (BBM2-4.  tAAtCC) )  / <2tAA> 

XX2-  (-BB-SQRT  (BBM2-4.  4AAICC) )  /  <2*AA) 

H0TX=0. 

IF(X(I+1)  .GE.  X ( I ) ) THEN 

CALL  LINSEG (XCI) , TFALL (I) , XXI, XX2, X (1+1) , TFALL ( 1+1 ) , HOTX) 
ELSE 

CALL  LINSEG <  X  < 1+1 ) , TFALL ( 1+ 1 ) , XX 1 , XX2, X  < I ) , TFALL <  I ) , HOTX ) 
END  IF 

IF (HOTX  .NE.  0. )THEN 
HOTY*A ( I ) IHOTX+B ( I ) 

GOTO  20 
ENDIF 

10  CONTINUE 

SWI».TRUE. 

RETURN 

20  CONTINUE 

T*(HOTX-X (I) ) / (X(I+1)-X(I) ) 

HOTTA-TtTFALL ( 1+1 ) +TFALL <1)1(1. -T> 

HOTSX-TtSHRX ( 1+1 ) +SHRX < I ) *  < 1 . -T) 

HOTSY«T*SHRY ( 1+1 ) +SHRY ( I >  *  < 1 . -T) 

END 

SUBROUTINE  SHEAR (X , Y, XX , YY, TA, DINF, UTRDR, DR1 , DR7, DR14, YM,  SO,  SZ# 

C  C , SHRX , SHRY , DR30 , NBURS , HLINE, WLINE, ALPHA  f 

C  BETA, FV , FF , CSW I , BLAT ) 

mmmmmtmmmtmtmmmtmmtmtmmmmmtm 

$  THIS  ROUTINE  CALCULATES  THE  DOSE  AT  (XX,YY>  FROH  THE  HOTLINE  (X,Y) 


VARIABLES 

X,Y  ■  HOTLINE  COORDINATES 

XX, YY  -  OFF  HOTLINE  COORDINATES 

TA  •  TIME  OF  ARRIVAL  IN  HOURS  AT  <X,Y> 

DINF  ■  THE  TIME  INTEGRATED  DOSE  TO  INFINITY  AT  (XX, YY) 
UTRDR  -  UNIT  TIME  REFERENCE  DOSE  RATE  AT  (XX, YY) 

DR1  -  THE  DOSE  RATE  AT  1  DAY  AT  <XX,YY) 

DR7  ■  THE  DOSE  RATE  AT  7  DAYS  AT  <XX,YY> 

DR14  -  THE  DOSE  RATE  AT  14  DAYS  AT  <XX,YY) 

DR30  »  THE  DOSE  RATE  AT  30  DAYS  AT  (XX, YY) 

YM  -  YIELD  IN  MEGATONS 

SO  -  INITIAL  STABILIZED  CLOUD  STANDARD  DEVIATION  IN  THE 
HORIZONTAL 


*  SZ  ■  STABILIZED  CLOUD  STANDARD  DEVIATION  IN  THE  VERTICAL 

*  C  -  COLARCO’S  COEFFICIENTS  FOR  R(T) 

*  SHRX  -  WIND  SHEAR  IN  THE  X  DIRECTION  AT  (X,Y> 

*  SHRY  -  WIND  SHEAR  IN  THE  Y  DIRECTION  AT  (X,Y> 

*  HC  «  STABILIZED  CLOUD  HEIGHT 

*  NBURS  -  NUMBER  OF  BURSTS 

«  HL I NE  ■  HEIGHT  OF  AREA  OF  BURSTS 

I  WLINE  -  WIDTH  OF  AREA  OF  BURSTS 

t  TASTAR, TC  -  TIME  VARIABLES  TO  FIND  SIGX,SIGY  ACCORDING  TO  WSEG 
«  SIGY  >  THE  STABILIZED  CLOUD  Y  STANDARD  DEVIATION  AT  <X,Y) 

*  SIGX  «  THE  STABILIZED  CLOUD  X  STANDARD  DEVIATION  AT  <X,Y) 

*  R  «  PARTICLE  RADIUS  IN  MICRONS  OR  METERS 

*  T  *  DUMMY  TIME  VARIABLE  TO  AVOID  GTA  SINGULARITY  AT  T«0 

*  RT  *  DR/DT  OF  FALLOUT  AT  (X,Y) 

*  ALPHA, BETA, ALPH2,ALPH3  *  PARAMETERS  FOR  ACTIVITY  SIZE 

*  DISTRIBUTION 

*  FV  «  ACTIVITY  VOLUME  FRACTION 
«  FF  «  FISSION  FRAVTION  OF  YIELD 

*  DENOM, CONST 1 , C0NST2  -  DUMMY  VARIABLES  TO  SIMPLIFY  EXPRESSIONS 

*  FOR  ACTIVITY  SIZE  DISTRIBUTION 

*  Z3,Z2  «  EXPONENTIAL  POWERS  FOR  ACTIVITY  SIZE  VOLUME  AND  AREA 

*  FRACTIONS 

*  AR  -  ACTIVITY  FOR  PARTICLES  OF  SIZE  R 

*  C1F  >  TERM  IN  HOPKINS’  SMEARING  EQUATION 

*  Z1  -  EXPONENTIAL  TERM  IN  HOPKINS’  SMEARING  EQUATION 

*  FXY  ■  THE  2  DIMENSIONAL  GAUSSIAN  TERM  OF  HOPKINS’ SMEARING 

*  EQUATION 

*  K  -  SOURCE  NORMALIZATION  CONSTANT 

*  GTA  -  THE  FRACTION  OF  ACTIVITY  ARRIVING  EVERYWHERE  AT  TA 

*  BLAT  »  BURST  LATITUDE 

*  CSWI  -  LOGICAL  SWITCH,  IF  .TRUE.  THEN  FINDING  DOSE  AT  A  CITY 


SUBROUTINES  CALLED 

SINGLE  -  FINDS  THE  AVERAGE  DOSE  ACROSS  A  CELL  OR  CITY  FOR  A 
SINGLE  BURST 

MULTI  -  FINDS  THE  DOSE  AT  A  CELL  CENTER  OR  CITY  FOR  AN  AREA  OF 
BURSTS 

NO  FILES  ARE  USED 

PARAMETER (PI-3. 14159265) 

REAL  X, Y,XX, YY,TA,DINF,UTRDR,DR1,DR7,DR14, YM,S0,SZ,C(7) ,SHRX, 

C  SHR Y , HC , S I GX , S I GY , R , T , RT , ALPHA , BET A , ALPH3 , ALPH2 , FV , DENOM , 

C  C0NST2, C0NST3, Z2,  Z3, AR, GTA, CIF, Zi, FXY, DR30,K,FF, NBURS, 

C  HLIN£,WLINE,D,B,A,Z6, Z7, ANG, ANGM 

LOGICAL  CSWI 
INTEGER  I 

DATA  K/7452. / 

HC-SZ/. 18*5. 280/ 1.609 
TASTAR»MIN(TA,3. ) 

TC-12. IHC/60. -2. 51 (HC/60. ) *42 


SIGY«S0tt2t  <1+8. tTASTAR/TC)  +  (SZtSHRYtTA/ 10. ) tt2 
SIGY-SQRT(SIGY) 

SIGX-S0tt2t  <1+8. tTASTAR/TC) ♦ <SZtSHRXITA/10. ) tt2 
SIGX-SQRT  <SIGX) 

T»MAX(TA,.l) 

R-0. 

DO  10  1-1,6 
R-R+C ( I ) ITtt ( 1-6) 

10  CONTINUE 

R-R+C <7> /SORT <T) 

R«Rtl.E06 

RT«0. 

DO  20  1-1,5 

RT-RT+ < 1-6) tC  < I ) tTt t  < 1-7) 

20  CONTINUE 

RT-RT-. 5IC (7) ITtt (-1.5) 

RT-RTtl.E06 

ALPH3-ALPHA+3 . IBETAtt2 
ALPH2-ALPHA+2 . t  BET At  1 2 
DENOM-SQRT  <2tPI) tBETAtR 
CONST3-FV/DENOM 
C0NST2- < 1 . -FV) /DENOM 

IF  <FV. GT. 1 . ) THEN 
ALPH3-ALPHA+2 . 5tBETAtt2 
FV-1. 

ENDIF 

Z3- <  <LOG <R) -ALPH3) /BETA) tt2 
Z2- <  <  LOG  <  R ) -ALPH2 ) /BET A ) 1 1 2 
AR-CQNST3tEXP  <-.  5tZ3>  +C0NST2tEXP  <-. 5tZ2> 

GTA— ARtRT 

IF<TA.LT.  . 1 )GTA=GTAtTA/. 1 
IF (NBURS  .LT.  2.) THEN 

CALL  SINGLE (FXY, X,Y, XX,  YY,SIGX,SIGY,BLAT,CSWI) 

ELSE 

CALL  MULTI <FXY, X, Y, XX, YY,SIGX, SI6Y, HLINE,MLINE, NBURS) 
ENDIF 


UTRDR-KtYMtFFtlOOO. tGTAtFXYtTA 
IF <UTRDR  .LT.  . lE-05)UTRDR-0. 
DINF-5. tUTRDRtTAtt <-. 2) 
DRl-UTRDRt24.tt  <-1.2) 
DR7-UTRDRt 168. tt  <-l . 2) 
DR14-UTRDRt336. tt (-1.2) 
DR30-UTRDRt720 . tt (-1 . 2) 

END 
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SUBROUTINE  SINGLE (FXY,X,Y, XX, YY,SI6X,SIGY, BLAT, CSWI) 
********************************************************************* 
FINDS  THE  AVERAGE  DOSE  ACROSS  A  CELL  OR  CITY  FROM  A  SINGLE  BURST 


VARIABLES 

FXY  -  THE  2D  GAUSSIAN  TERM  IN  HOPKINS’  SMEARING  EQUATION 

X,Y  ■  HOTLINE  COORDINATE 

XX, YY  *  COORDINATE  OF  A  CITY  OR  CELL  CENTER 

SIGX,SIGY  -  STANDARD  DEVIATION  OF  2D  GAUSSIAN  IN  X,Y  DIRECTION 

BLAT  *  BURST  LATITUDE 

CSWI  -  LOGICAL  SWITCH,  IF  .TRUE.  THEN  FINDING  DOSE  AT  A  CITY 
Z1  -  DISTANCE  IN  STANDARD  DEVIATIONS  SQUARED  FROM  (X,Y>  TO 
(XX, YY) 

D  -  DISTANCE  FROM  <X,Y>  TO  (XX, YY) 

HLINE, WLINE  ■  HEIGHT  AND  WIDTH  OF  THE  ONE  DEGREE  BY  ONE  DEGREE 
LATITUDE  BY  LONGTITUDE  CELL 
ANGM  -  ANGLE  FROM  THE  Y  AXIS  TO  THE  CELL  DIAGONAL 
ANG  -  ANGLE  FROM  THE  X  AXIS  TO  THE  EFFECTIVE  WIND  VECTOR 
A,B  -  LIMITS  OF  INTEGRATION  TO  FIND  AVERAGE  DOSE  ACROSS  A  CELL 
OR  CITY,  IN  KM 

UL,LL  =  LIMITS  OF  INTEGRATION  TO  FIND  AVERAGE  DOSE,  IN  STANDARD 
DEVIATIONS 


SUBROUTINES  CALLED 

CNF  -  CUMULATIVE  NORMAL  FUNCTION  TO  FIND  AVERAGE  DOSE 


NO  FILES  ARE  USED 

********************************************************************* 


PARAMETER (PI-3. 14159265) 

REAL  FXY,  X, Y, XX, YY, SIGX, SIGY, BLAT, Z 1 , D, HLINE, WLINE, ANGM, ANG, B, 
C  A,SIG,UL,LL 
LOGICAL  CSWI 

Zl-(XXIY-YYtX) $12/ ( (SIGYtX) **2+(SIGX*Y) 1*2) 

IF (Z1  .GT.  500. )THEN 
FXY-O. 

RETURN 

ENDIF 

D-SQRT ( (XX-X) 1*2+ (YY-Y) **2) 

HLINE-110.94 

WLINE-C0S(BLAT*PI/180. ) $110. 94 
ANGM-ATAN (WLINE/HLINE) 

ANG-ATAN(Y/X) 

IF (ABS(ANG)  .GT.  ANGM) THEN 
B-D+. 5IWLINE/ABS (SIN (ANG) ) 

A-D-.5IWLINE/AB8 (SIN (ANG) ) 

ELSE 

B-D+ . 5*HL I NE /COS ( ANG ) 

A-D-. 5IHLINE/C0S (ANG) 


IP (CSWl)THEN 
B“D+38. 

A-D-38. 

END  IF 

SIG-SQRT( (SIGYIX) **2+(SIGX*Y) t*2) 

UL*  < -B*Y$SIN ( ANG) -B*X*COS (ANG) ) /SIG 
LL« (-AIYISIN (ANG) -AlXtCOS ( ANG) ) /SIG 
FXY- (CNF (UL) -CNF (LL) > / ( (B-A) * (-Y*SIN (ANG) -X*CQS (ANG) > ) 

END 

SUBROUTINE  MULTI (FXY, X, Y, XX, YY, SIGX, SIGY, HLINE, WLINE, NBURS) 

ttttttttttittttttttttttttitttttttttttttttiititttttttttttttttttttttttt 

i  THIS  ROUTINE  FINDS  THE  DOSE  AT  A  CELL  CENTER  OR  CITY  FOR  AN  AREA 
OF  BURSTS 


VARIABLES 

FXY  *  2D  GAUSSIAN  TERM  IN  HOPKINS*  SMEARING  EQUATION 

X,Y  «  HOTLINE  COORDINATE 

XX, YY  =  CELL  CENTER  OR  CITY  COORDINATE 

SIGX,SIGY  *  STANDARD  DEVIATION  OF  2D  GAUSSIAN  CLOUD  IN  X,Y 
DIRECTION 

HLINE, WLINE  *  HEIGHT  AND  WIDTH  OF  AN  AREA  OF  BURSTS 
NBURS  ■=  NUMBER  OF  BURSTS 

Z1  -  DISTANCE  IN  STANDARD  DEVIATIONS  FROM  (X,Y>  TO  (XX, YY) 

ANGM  «  ANGLE  FROM  X  AXIS  TO  DIAGONAL  OF  AREA  OF  BURSTS 
ANG  =  ANGLE  BETWEEN  X  AXIS  AND  EFFECTIVE  WIND  VECTOR 
DIA  =  DIAGONAL  OF  THE  AREA  OF  BURSTS 
L  -  BASE  OF  BURST  DENSITY  TRAPEZOID 
B  ■  TOP  OF  BURST  DENSITY  TRAPEZOID 
H  -  HEIGHT  OF  BURST  DENSITY  TRAPEZOID 

SIGP  *  KM  PER  STANDARD  DEVIATION  IN  THE  CROSSWIND  DIRECTION 
SIGL  -  NUMBER  OF  STANDARD  DEVIATIONS  IN  CROSSWIND  DIRECTION 
FROM  CENTER  OF  TRAPEZOID  BASE  TO  END  OF  BASE 
FACT  *  CONSTANT  FACTOR  FROM  LINE  INTEGRATION  OF  TRAPEZOID 
SIG  -  CONSTANT  FACTOR  FROM  LINE  INTEGRATION  OF  TRAPEZOID 
UL,LL  -  LIMITS  OF  INTEGRATION  IN  NUMBER  OF  STANDARD  DEVIATIONS 
UL1.LL1  -  DUMMY  VARIABLES  TO  PREVENT  TAKING  EXPONENT  OF  TOO 
SMALL  OF  A  NUMBER 

MULTI, MULT2,MULT3  -  CONSTANTS  IN  FRONT  OF  INTEGRALS  FOR  ENDS  OF 

TRAPEZOID  INTEGRATION 

PARTI , PART2, PART3  -  THREE  PARTS  OF  EACH  INTEGRAL  ON  THE  ENDS  OF 

THE  TRAPEZOID 

INTI, INT2, INT3  «  INTEGRALS  FOR  THE  THREE  REGIONS  OF  THE 
TRAPEZOID 


SUBROUTINES  CALLED 

CNF  -  CUMULATIVE  NORMAL  FUNCTION 


*  NO  FILES  ARE  USED 


PARAMETER (PI-3. 14139263) 

REAL  FXY, X,Y, XX, YY,SIGX,SIGY,HLINE, WLINE,NBURS, Zl, ANGM, AN6,DIA 
L,B,H,SIGP,SIGL, FACT, SIG, UL,LL, MULTI, MULT2,MULT3, PARTI, 

I  PART2,PART3, INTI, INT2,  INT3,C0T 
COT (X)-1./TAN(X) 

Zl-SQRT ( (XXtY-YYtX) 1*2/ ( (SIGYtX) tt2+ (SIGXtY) tt2> > 

ANGM® AT  AN ( HL I NE / WL I NE ) 

ANG-ABS(ATAN(Y/X) ) 

DIA=SQRT(HLINEtt2+WLINEtt2) 

IF(AN6  .LE.  ANGM) THEN 
L=DIAtSIN (ANG+ANGM) 

B=COS(ANG) t (HLINE-WLINEtTAN (ANG) ) 

H-2. tNBURS/ (B+L) 

ELSE 

L=DIAtSIN (ANG+ANGM) 

B-SIN(ANG) « ( WL I NE-HL I  NE  t  COT (ANG) > 

H-2. tNBURS/ (B+L) 

ENDIF 

ANG=ATAN(Y/X) 

SIGP=SQRT ( (SIGXtSIN(ANG) ) 4>2+ (SIGYtCOS(ANG) ) %%2) 

SIGL-L/ (2. tSIGP) 

IF(Z1  .GT.  (SIGL+3. ) )THEN 
FXY=0. 

RETURN 

ENDIF 

FACT-YtSIN ( ANG) +X*COS (ANG) 

SIG-SQRT ( (SIGXtY) It2+ (SI GY*  X ) tt2) 

IF ( (L-B)  .LT.  L/IOO. )THEN 
L-(L+B)/2. 

UL- ( ( XX+L/2. tSIN (ANG) ) tY- (YY-L/2. tCOS (ANG) ) IX) /SIG 
LL- ( (XX-L/2. ISIN (ANG) ) IY- ( YY+L/2. tCOS (ANG) ) tX) /SIG 
FXY-NBURS/ (LtFACT) t (CNF (UL) -CNF (LL) ) 

RETURN 

ENDIF 

UL«( (XX-B/2. tSIN(ANG) ) tY-(YY+B/2. tCOS (ANG) ) tX) /SIG 
LL- ( (XX-L/2. tSIN (ANG) >tY- (YY+L/2. tCOS (ANG) )tX) /SIG 
MULTI— 2. tHtSIG/ ( (L-B) tFACTtt2tSQRT  <2. tPI ) ) 

UL1-MIN (ABS (UL) ,30. ) 

LL1-MIN(ABS (LL) , 30. ) 

PARTI -MULT It (EXP (-.5tULltt2) -EXP(-.5tLLltt2) > 

MULT2-2. tHt (YYtX-XXtY) / ( (L-B) tFACTtt2) 

PART2-MULT2t (CNF (UL)-CNF (LL) ) 

MULT3-LtH/ ( (L-B) tFACT) 

PART3-MULT3t (CNF (UL) -CNF (LL) ) 

I NT 1 -PART 1 +PART2+PART3 


UL« ( ( XX+B/2. tSIN (ANG) ) IY- < YY-B/2. ICOS < ANG) ) I X  > /SIG 
LL* < (XX-B/2. tSIN (ANG) ) tY- (YY+B/2. ICOS (ANG) ) IX) /SIG 
INT2*H/FACTt (CNF (UL) -CNF (LL) ) 

UL* ( <  XX+L/2. tSIN (ANG) ) IY- < YY-L/2. ICOS (ANG) >  IX) /SIG 
LL= < (XX+B/2. tSIN (ANG) ) tY- (YY-B/2. ICOS (ANG) ) *X) /SIG 
UL1=MIN(ABS(UL) ,30. ) 

LL1*MIN(ABS(LL),30.) 

PART1=-MULT1 t (EXP (-. 5IUL1 112) -EXP (-. 5ILL1 112) ) 
PART2=-MULT2* (CNF (UL) -CNF (LL) ) 

PART3«MULT3* (CNF (UL) -CNF (LL) ) 

I NT3*PART 1 +PART2+PART3 

FXY=INT1+INT2+INT3 

END 


SUBROUTINE  SLQINT (X, Y, A, B, NP) 

tttttttttttttuttttuttttttttttttttttttttuttttttutttttttttttttttttt 

THIS  ROUTINE  FINDS  THE  SLOPE  AND  INTERCEPT  OF  EACH  LINE  SEGMENT 
DEFINED  BY  THE  HOTLINE 


VARIABLES 

X  «  THE  HOTLINE  X  COORDINATES  IN  KM 

Y  -  THE  HOTLINE  Y  COORDINATES  IN  KM 

A(I)  -  THE  SLOPE  OF  THE  HOTLINE  SEGMENT  FROM  I  TO  1+1 

B ( I )  «  THE  INTERCEPT  OF  THE  HOTLINE  SEGMENT 


*  NO  SUBROUTINES  ARE  CALLED 

$ - 

*  NO  FILES  ARE  USED 

REAL  X (0:NP) , Y (Ot NP) , A(Oi NP-1) ,B(OsNP-l) 

INTEGER  I 

DO  10  I-O.NP-l 

A(I)*(Y(I+1)-Y(I))/(X(I+1)-X(I)) 

B (I)*Y ( I)-A< I) IX (I) 

10  CONTINUE 

END 

SUBROUTINE  LINSEG ( XL ,  TAL , X X 1 , X X2 , XG , TAG, HOTX ) 

tmmimmttmitmmmtmmmmttmmtmmimtmtm 

*  THIS  ROUTINE  FINDS  WHETHER  XXI  AND/OR  XX2  IS  ON  THE  HOTLINE 
I  SEGMENT  DEFINED  BY  XL  AND  XG.  HOTX  IS  THE  CLOSEST  POINT  TO  THE 

I  ORIGIN  ON  THE  SEGMENT.  IF  NEITHER  POINT  IS  ON  THE  SEGMENT,  H0TX«0. 


VARIABLES 

XL  -  X  LESSER  OF  THE  TWO  ENDPOINTS  DEFINING  THE  SEGMENT 
XXI  -  ONE  OF  THE  QUADRATIC  ROOTS  ON  THE  LINE 


XX2  =  ONE  OF  THE  QUADRATIC  ROOTS  ON  THE  LINE 

XG  =  X  GREATER  OF  THE  TWO  ENDPOINTS  DEFINING  THE  LINE  SEGMENT 

HOTX  =  THE  POINT  XXI  OR  XX2  WHICH  IS  ON  THE  LINE  SEGMENT 

TAL  *  TIME  OF  ARRIVAL  FOR  XL,  THE  SMALLEST  X  ON  HOTLINE  SEGMENT 

TAG  =  TIME  OF  ARRIVAL  FOR  XG,  THE  GREATEST  X  ON  HOTLINE  SEGMENT 

NO  SUBROUTINES  ARE  CALLED 

*  NO  FILES  ARE  USED 

REAL  XL, XXI, XX2,XG, HOTX, TAL, TAG 
IFtXXl  .EQ.  0. ) X X 1 = 1 . E25 
IF (XX2  .EQ.  0.)XX2=1.E25 

IF (XXI  .GE.  XL  .AND.  XXI  .LE.  XG)H0TX=XX1 
IF ( XX2  .GE.  XL  .AND.  XX2  .LE.  XG)H0TX=XX2 

IF (XXI  .GE. XL  .AND.  XX2  .GE.  XL  .AND.  XXI  .LE.  XG  .AND. 

C  XX2  .LE.  XG)THEN 

IF (TAL  .LT.  TAG) THEN 
IF(ABS(XX1-XL)  .LT.  ABS(XX2-XL; >THEN 
H0TX*XX1 
ELSE 

H0TX=XX2 
END  IF 

ELSE 

IF (ABS(XXl-XG)  .LT.  ABS(XX2-XG) )THEN 
H0TX=XX1 
ELSE 

H0TX=XX2 
ENDIF 
END  IF 
ENDIF 

END 

REAL  FUNCTION  CNF(Z) 

*  CALCULATES  THE  CUMULATIVE  NORMAL  FUNCTION,  BY  POLYNOMIAL  APPROX. 

I  FROM  ABRAMOWITZ  AND  STEGUN. 


VARIABLES 

C1,C2,C3,C4  =  POLYNOMIAL  COEFFICIENTS 
Z  «  ARGUMENT  OF  FUNCTION 
X  -  VARIABLE  IN  POLYNOMIAL  EXPANSION 


NO  SUBROUTINES  ARE  CALLED 


*  NO  FILES  ARE  USED 


1.0  X 


Itf  1£8  121 


Hi  |H  Boo 

iu  K±± 

iu  |££ 

!  ^  Ih 


L25  IU  1 1.6 


MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BUREAU  OF  STANDARDS- 1963-A 


PARAMETER  <C1«. 196854, C 2- . 1 15194, C3-. 000344, C4-. 019259) 
REAL  X,Z 

X-ABS(Z) 

CNF-1 . +C1 »X+C2*Xtf2+C3*X«*3+C4*X*»4 

CNF-. 5/CNF* *4 

IF ( Z. BE. 0.) CNF-1. -CNF 
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‘Hopkins'  variable  wind  fallout  model  is  used  to  predict  the  dose  and 
population  insult  across  the  United  States  from  a  nuclear  attack.  The 
dose  calculation  is  performed  by  two  programs  written  in  Fortran  V  for  a 
CYBER  845  computer.  Hopkins'  hotline  locator  program  was  modified  to 
reduce  its  run  time,  and  it  is  used  to  locate  the  fallout  hotline  as 
trace  particles  are  translated  to  the  ground  in  a  spatially  varying  wind 
field.  The  second  program  analytically  smears  fallout  activity  along  the 
hotline.  To  reduce  run  time  and  to  match  the  population  model,  the  dose 
program  uses  a  computational  grid  of  one  degree  latitude  by  one  degree 
longitude.  A  difference  of  cumulative  normal  functions  gives  the  average 
dose  across  a  grid  cell.  An  analytical  method  was  developed  to  treat 
multiple  bursts  against  an  area  target  as  one  cloud. 

For  the  winds  of  0000  Universal  Time  on  16  January  1982,  a  hypo¬ 
thetical  attack  against  twenty  five  air  bases  and  six  Minuteman  missile 
fields  results  in  26.9  million  fallout  deaths.  This  calculation  used 
407  seconds  of  computer  time.  2 
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